The R Project SVN – R – Diff – /branches/djm-source/src/main/plot3d.c


/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 *  Copyright (C) 1997-2002   Robert Gentleman, Ross Ihaka and the R core team.
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <Defn.h>
#include <Rmath.h>
#include <Graphics.h>
#include <Rdevices.h>
#include <Print.h>
#include <R_ext/Boolean.h>
 
#ifndef HAVE_HYPOT
# define hypot pythag
#endif
 
/* Conversion of degrees to radians */
 
#define DegToRad(x) (0.01745329251994329576 * x)
 
/* Definitions of data structures for vectors and */
/* transformations in homogeneous 3d coordinates */
 
typedef double Vector3d[4];
typedef double Trans3d[4][4];
 
/* The viewing transformation matrix. */
 
static SEXP gcall;
static Trans3d VT;
 
static void TransVector (Vector3d u, Trans3d T, Vector3d v)
{
    double sum;
    int i, j;
 
    for (i = 0; i < 4; i++) {
	sum = 0;
	for (j = 0; j < 4; j++)
	    sum = sum + u[j] * T[j][i];
	v[i] = sum;
    }
}
 
static void Accumulate (Trans3d T)
{
    Trans3d U;
    double sum;
    int i, j, k;
 
    for (i = 0; i < 4; i++) {
	for (j = 0; j < 4; j++) {
	    sum = 0;
	    for (k = 0; k < 4; k++)
		sum = sum + VT[i][k] * T[k][j];
	    U[i][j] = sum;
	}
    }
    for (i = 0; i < 4; i++)
	for (j = 0; j < 4; j++)
	    VT[i][j] = U[i][j];
}
 
static void SetToIdentity (Trans3d T)
{
    int i, j;
    for (i = 0; i < 4; i++) {
	for (j = 0; j < 4; j++)
	    T[i][j] = 0;
	T[i][i] = 1;
    }
}
 
static void Translate (double x, double y, double z)
{
    Trans3d T;
    SetToIdentity(T);
    T[3][0] = x;
    T[3][1] = y;
    T[3][2] = z;
    Accumulate(T);
}
 
static void Scale (double x, double y, double z)
{
    Trans3d T;
    SetToIdentity(T);
    T[0][0] = x;
    T[1][1] = y;
    T[2][2] = z;
    Accumulate(T);
}
 
static void XRotate (double angle)
{
    double c, s;
    Trans3d T;
    SetToIdentity(T);
    c = cos(DegToRad(angle));
    s = sin(DegToRad(angle));
    T[1][1] = c;
    T[2][1] = -s;
    T[2][2] = c;
    T[1][2] = s;
    Accumulate(T);
}
 
static void YRotate (double angle)
{
    double c, s;
    Trans3d T;
    SetToIdentity(T);
    c = cos(DegToRad(angle));
    s = sin(DegToRad(angle));
    T[0][0] = c;
    T[2][0] = s;
    T[2][2] = c;
    T[0][2] = -s;
    Accumulate(T);
}
 
static void ZRotate (double angle)
{
    double c, s;
    Trans3d T;
    SetToIdentity(T);
    c = cos(DegToRad(angle));
    s = sin(DegToRad(angle));
    T[0][0] = c;
    T[1][0] = -s;
    T[1][1] = c;
    T[0][1] = s;
    Accumulate(T);
}
 
static void Perspective (double d)
{
    Trans3d T;
 
    SetToIdentity(T);
    T[2][3] = -1 / d;
    Accumulate(T);
}
 
 
/* Stuff for labels on contour plots
   Originally written by Nicholas Hildreth
   Adapted by Paul Murrell
*/
 
static
void FindCorners(double width, double height, SEXP label,
		 double x0, double y0, double x1, double y1,
		 DevDesc *dd) {
    double delta = height / width;
    double dx = GConvertXUnits(x1 - x0, USER, INCHES, dd) * delta;
    double dy = GConvertYUnits(y1 - y0, USER, INCHES, dd) * delta;
    dx = GConvertYUnits(dx, INCHES, USER, dd);
    dy = GConvertXUnits(dy, INCHES, USER, dd);
 
    REAL(label)[0] = x0 + dy;
    REAL(label)[4] = y0 - dx;
    REAL(label)[1] = x0 - dy;
    REAL(label)[5] = y0 + dx;
    REAL(label)[3] = x1 + dy;
    REAL(label)[7] = y1 - dx;
    REAL(label)[2] = x1 - dy;
    REAL(label)[6] = y1 + dx;
}
 
static
int TestLabelIntersection(SEXP label1, SEXP label2) {
 
    int i, j, l1, l2;
    double Ax, Bx, Ay, By, ax, ay, bx, by, q1, q2;
    double dom;
    double result1, result2;
 
    for (i = 0; i < 4; i++) {
	Ax = REAL(label1)[i];
	Ay = REAL(label1)[i+4];
	Bx = REAL(label1)[(i+1)%4];
	By = REAL(label1)[(i+1)%4+4];
	for (j = 0; j < 4; j++) {
	    ax = REAL(label2)[j];
	    ay = REAL(label2)[j+4];
	    bx = REAL(label2)[(j+1)%4];
	    by = REAL(label2)[(j+1)%4+4];
 
	    q1 = Ax*(ay-by);
	    q2 = Ay*(bx-ax);
 
	    dom = Bx*by - Bx*ay - Ax*by + Ax*ay - bx*By + bx*Ay + ax*By - ax*Ay;
	    if (dom == 0.0) {
		result1 = -1;
		result2 = -1;
	    }
	    else {
		result1 = (bx*Ay - ax*Ay - ay*bx - Ax*by + Ax*ay + by*ax) / dom;
 
		if (bx - ax == 0.0) {
		    if (by - ay == 0.0)
			result2 = -1;
		    else
			result2 = (Ay + (By - Ay) * result1 - ay) / (by - ay);
		}
		else
		    result2 = (Ax + (Bx - Ax) * result1 - ax) / (bx - ax);
 
	    }
	    l1 = (result1 >= 0.0) && (result1 <= 1.0);
	    l2 = (result2 >= 0.0) && (result2 <= 1.0);
	    if (l1 && l2) return 1;
	}
    }
 
    return 0;
}
 
/*** Checks whether a label window is inside view region ***/
static int LabelInsideWindow(SEXP label, DevDesc *dd) {
    int i = 0;
    double x, y;
 
    while (i < 4) {
	x = REAL(label)[i];
	y = REAL(label)[i+4];
	GConvert(&x, &y, USER, NDC, dd);
	/*	x = GConvertXUnits(REAL(label)[i], USER, NDC, dd);
		y = GConvertYUnits(REAL(label)[i+4], USER, NDC, dd); */
 
	if ((x < 0) || (x > 1) ||
	    (y < 0) || (y > 1))
	    return 1;
	i += 1;
    }
    return 0;
}
 
 
	/*  C o n t o u r   P l o t t i n g  */
 
typedef struct SEG {
    struct SEG *next;
    double x0;
    double y0;
    double x1;
    double y1;
} SEG, *SEGP;
 
static SEGP *ctr_SegDB;
 
static int ctr_intersect(double z0, double z1, double zc, double *f)
{
    if ((z0 - zc) * (z1 - zc) < 0.0) {
	*f = (zc - z0) / (z1 -	z0);
	return 1;
    }
    return 0;
}
 
static SEGP ctr_newseg(double x0, double y0, double x1, double y1, SEGP prev)
{
    SEGP seg = (SEGP)R_alloc(1, sizeof(SEG));
    seg->x0 = x0;
    seg->y0 = y0;
    seg->x1 = x1;
    seg->y1 = y1;
    seg->next = prev;
    return seg;
}
 
static void ctr_swapseg(SEGP seg)
{
    double x, y;
    x = seg->x0;
    y = seg->y0;
    seg->x0 = seg->x1;
    seg->y0 = seg->y1;
    seg->x1 = x;
    seg->y1 = y;
}
 
	/* ctr_segdir(): Determine the entry direction to the next cell */
	/* and update the cell indices */
 
#define XMATCH(x0,x1) (fabs(x0-x1) == 0)
#define YMATCH(y0,y1) (fabs(y0-y1) == 0)
 
static int ctr_segdir(double xend, double yend, double *x, double *y,
		      int *i, int *j, int nx, int ny)
{
    if (YMATCH(yend, y[*j])) {
	if (*j == 0)
	    return 0;
	*j = *j - 1;
	return 3;
    }
    if (XMATCH(xend, x[*i])) {
	if (*i == 0)
	    return 0;
	*i = *i - 1;
	return 4;
    }
    if (YMATCH(yend, y[*j + 1])) {
	if (*j >= ny - 1)
	    return 0;
	*j = *j + 1;
	return 1;
    }
    if (XMATCH(xend, x[*i + 1])) {
	if (*i >= nx - 1)
	    return 0;
	*i = *i + 1;
	return 2;
    }
    return 0;
}
 
/* Search seglist for a segment with endpoint (xend, yend). */
/* The cell entry direction is dir, and if tail=1/0 we are */
/* building the tail/head of a contour.	 The matching segment */
/* is pointed to by seg and the updated segment list (with */
/* the matched segment stripped) is returned by the funtion. */
 
static SEGP ctr_segupdate(double xend, double yend, int dir, Rboolean tail,
			  SEGP seglist, SEGP* seg)
{
    if (seglist == NULL) {
	*seg = NULL;
	return NULL;
    }
    switch (dir) {
    case 1:
    case 3:
	if (YMATCH(yend,seglist->y0)) {
	    if (!tail)
		ctr_swapseg(seglist);
	    *seg = seglist;
	    return seglist->next;
	}
	if (YMATCH(yend,seglist->y1)) {
	    if (tail)
		ctr_swapseg(seglist);
	    *seg = seglist;
	    return seglist->next;
	}
	break;
    case 2:
    case 4:
	if (XMATCH(xend,seglist->x0)) {
	    if (!tail)
		ctr_swapseg(seglist);
	    *seg = seglist;
	    return seglist->next;
	}
	if (XMATCH(xend,seglist->x1)) {
	    if (tail)
		ctr_swapseg(seglist);
	    *seg = seglist;
	    return seglist->next;
	}
	break;
    }
    seglist->next = ctr_segupdate(xend, yend, dir, tail, seglist->next, seg);
    return seglist;
}
 
/* labelList, label1, and label2 are all SEXPs rather than being allocated
   using R_alloc because they need to persist across calls to contour().
   In do_contour() there is a vmaxget() ... vmaxset() around each call to
   contour() to release all of the memory used in the drawing of the
   contour _lines_ at each contour level.  We need to keep track of the
   contour _labels_ for _all_ contour levels, hence we have to use a
   different memory allocation mechanism.
*/
static SEXP labelList;
 
static
double distFromEdge(double *xxx, double *yyy, int iii, DevDesc *dd) {
    return fmin2(fmin2(xxx[iii]-Rf_gpptr(dd)->usr[0], Rf_gpptr(dd)->usr[1]-xxx[iii]),
		 fmin2(yyy[iii]-Rf_gpptr(dd)->usr[2], Rf_gpptr(dd)->usr[3]-yyy[iii]));
}
 
static
Rboolean useStart(double *xxx, double *yyy, int ns, DevDesc *dd) {
    if (distFromEdge(xxx, yyy, 0, dd) < distFromEdge(xxx, yyy, ns-1, dd))
	return TRUE;
    else
	return FALSE;
}
 
static
int findGapUp(double *xxx, double *yyy, int ns, double labelDistance,
	      DevDesc *dd) {
    double dX, dY;
    double dXC, dYC;
    double distanceSum = 0;
    int n = 0;
    int jjj = 1;
    while ((jjj < ns) && (distanceSum < labelDistance)) {
	/* Find a gap big enough for the label
	   use several segments if necessary
	*/
	dX = xxx[jjj] - xxx[jjj - n - 1]; /* jjj - n - 1 == 0 */
	dY = yyy[jjj] - yyy[jjj - n - 1];
	dXC = GConvertXUnits(dX, USER, INCHES, dd);
	dYC = GConvertYUnits(dY, USER, INCHES, dd);
	distanceSum = hypot(dXC, dYC);
	jjj++;
	n++;
    }
    if (distanceSum < labelDistance)
	return 0;
    else
	return n;
}
 
static
int findGapDown(double *xxx, double *yyy, int ns, double labelDistance,
		DevDesc *dd) {
    double dX, dY;
    double dXC, dYC;
    double distanceSum = 0;
    int n = 0;
    int jjj = ns - 2;
    while ((jjj > -1) && (distanceSum < labelDistance)) {
	/* Find a gap big enough for the label
	   use several segments if necessary
	*/
	dX = xxx[jjj] - xxx[jjj + n + 1]; /*jjj + n + 1 == ns -1 */
	dY = yyy[jjj] - yyy[jjj + n + 1];
	dXC = GConvertXUnits(dX, USER, INCHES, dd);
	dYC = GConvertYUnits(dY, USER, INCHES, dd);
	distanceSum = hypot(dXC, dYC);
	jjj--;
	n++;
    }
    if (distanceSum < labelDistance)
	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, names;
    /* 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);
		/* 
		 * Set the names attribute for the contour
		 * So that users can extract components using 
		 * meaningful names
		 */
		PROTECT(names = allocVector(STRSXP, 3));
		SET_STRING_ELT(names, 0, mkChar("level"));
		SET_STRING_ELT(names, 1, mkChar("x"));
		SET_STRING_ELT(names, 2, mkChar("y"));
		setAttrib(ctr, R_NamesSymbol, names);
		/* 
		 * 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(5);
	    }
	}
    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();
 
SEXP do_contourLines(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP oargs, c, x, y, z;
    int nx, ny, nc;
    GEDevDesc *dd = GEcurrentDevice();
 
    oargs = args;
 
    x = CAR(args);
    internalTypeCheck(call, x, REALSXP);
    nx = LENGTH(x);
    args = CDR(args);
 
    y = CAR(args);
    internalTypeCheck(call, y, REALSXP);
    ny = LENGTH(y);
    args = CDR(args);
 
    z = CAR(args);
    internalTypeCheck(call, z, REALSXP);
    args = CDR(args);
 
    /* levels */
    c = CAR(args);
    internalTypeCheck(call, c, REALSXP);
    nc = LENGTH(c);
    args = CDR(args);
    
    return GEcontourLines(REAL(x), nx, REAL(y), ny, REAL(z), REAL(c), nc, dd);
}
 
/*
 * The *base graphics* function contour() and the *general base* 
 * function contourLines() use the same code to generate contour lines 
 * (i.e., the function contourLines())
 *
 * I had a look at extracting the code that draws the labels 
 * into a *general base* function
 * (e.g., into some sort of labelLines() function), 
 * but the code is too base-graphics-specific (e.g., one of the 
 * labelling methods seeks the location closest to the edge of the
 * plotting region) so I've left it alone for now. 
 * 
 * This does mean that the contourLines() function is part of the
 * graphics engine, but the contour() function is part of the 
 * base graphics system.
 */
 
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' */
 
    char *vmax;
 
    double xend, yend;
    int i, ii, j, jj, ns, ns2, dir;
    SEGP seglist, seg, s, start, end;
    double *xxx, *yyy;
 
    double variance, dX, dY, deltaX, deltaY;
    double dXC, dYC, deltaXC, deltaYC;
    int range=0, indx=0, n; /* -Wall */
    double lowestVariance;
    double squareSum, sum;
    int iii, jjj;
    double distanceSum, labelDistance, avgGradient;
    int zeroCount;
    char buffer[255];
    double avg;
    int result;
    double ux, uy, vx, vy;
    double xStart, yStart;
    double dx, dy, dxy;
    double labelHeight;
    SEXP label1 = PROTECT(allocVector(REALSXP, 8));
    SEXP label2;
    SEXP lab;
    Rboolean gotLabel = FALSE;
    Rboolean ddl;/* Don't draw label -- currently unused, i.e. always FALSE*/
 
#ifdef DEBUG_contour
    Rprintf("contour(lev = %g):\n", zc);
#endif
 
    vmax = vmaxget();
    ctr_SegDB = contourLines(REAL(x), nx, REAL(y), ny, REAL(z), zc, atom);
    /* we need to keep ctr_SegDB available, so vmaxset(vmax); was wrong */
 
    /* The segment database is now assembled. */
    /* Begin following contours. */
    /* 1. Grab a segment */
    /* 2. Follow its tail */
    /* 3. Follow its head */
    /* 4. Draw the contour */
 
    for (i = 0; i < nx - 1; i++)
      for (j = 0; j < ny - 1; j++) {
	while ((seglist = ctr_SegDB[i + j * nx])) {
	    ii = i; jj = j;
	    start = end = seglist;
	    ctr_SegDB[i + j * nx] = seglist->next;
	    xend = seglist->x1;
	    yend = seglist->y1;
	    while ((dir = ctr_segdir(xend, yend, REAL(x), REAL(y),
				     &ii, &jj, nx, ny))) {
		ctr_SegDB[ii + jj * nx]
		    = ctr_segupdate(xend, yend, dir, TRUE,/* = tail */
				    ctr_SegDB[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, REAL(x), REAL(y),
				     &ii, &jj, nx, ny))) {
		ctr_SegDB[ii + jj * nx]
		    = ctr_segupdate(xend, yend, dir, FALSE,/* ie. head */
				    ctr_SegDB[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;
 
	    vmax = vmaxget();
	    xxx = (double *) R_alloc(ns + 1, sizeof(double));
	    yyy = (double *) R_alloc(ns + 1, sizeof(double));
	    /* now have the space, go through again: */
	    s = start;
	    ns = 0;
	    xxx[ns] = s->x0;
	    yyy[ns++] = s->y0;
	    while (s->next && ns < MAX_ns) {
		s = s->next;
		xxx[ns] = s->x0;
		yyy[ns++] = s->y0;
	    }
	    xxx[ns] = s->x1;
	    yyy[ns++] = s->y1;
#ifdef DEBUG_contour
	    Rprintf("  [%2d,%2d]: (x,y)[1:%d] = ", i,j, ns);
	    if(ns >= 5)
		Rprintf(" (%g,%g), (%g,%g), ..., (%g,%g)\n", 
			xxx[0],yyy[0], xxx[1],yyy[1], xxx[ns-1],yyy[ns-1]);
	    else
		for(iii = 0; iii < ns; iii++)
		    Rprintf(" (%g,%g)%s", xxx[iii],yyy[iii],
			    (iii < ns-1) ? "," : "\n");
#endif	    
 
	    GMode(1, dd);
 
	    if (drawLabels) {
		/* If user supplied labels, use i'th one of them
		   Otherwise stringify the z-value of the contour */
		buffer[0] = ' ';
		if (!isNull(labels)) {
		    int numl = length(labels);
		    strcpy(&buffer[1], CHAR(STRING_ELT(labels, cnum % numl)));
		}
		else {
		    PROTECT(lab = allocVector(REALSXP, 1));
		    REAL(lab)[0] = zc;
		    lab = labelformat(lab);
		    strcpy(&buffer[1], CHAR(STRING_ELT(lab, 0)));
		    UNPROTECT(1);
		}
		buffer[strlen(buffer)+1] = '\0';
		buffer[strlen(buffer)] = ' ';
 
		if (vectorFonts) {
		    /* 1, 1 => sans serif, basic font */
		    labelDistance = 
			GVStrWidth((unsigned char *)buffer, typeface, fontindex,
				   INCHES, dd);
		    labelHeight = 
			GVStrHeight((unsigned char *)buffer, typeface, fontindex,
				    INCHES, dd);
		}
		else {
		    labelDistance = GStrWidth(buffer, INCHES, dd);
		    labelHeight = GStrHeight(buffer, INCHES, dd);
		}
 
		if (labelDistance > 0) {
		    /* Try to find somewhere to draw the label */
		    switch (method) {
		    case 0: /* draw label at one end of contour
			       overwriting contour line
			    */
			if (useStart(xxx, yyy, ns, dd) )
			    indx = 0;
			else
			    indx = ns - 1;
			break;
		    case 1: /* draw label at one end of contour
			       embedded in contour
			       no overlapping labels
			    */
			indx = 0;
			range = 0;
			gotLabel = FALSE;
			if (useStart(xxx, yyy, ns, dd)) {
			    iii = 0;
			    n = findGapUp(xxx, yyy, ns, labelDistance, dd);
			}
			else {
			    n = findGapDown(xxx, yyy, ns, labelDistance, dd);
			    iii = ns - n - 1;
			}
			if (n > 0) {
			    /** Find 4 corners of label extents **/
			    FindCorners(labelDistance, labelHeight, label1,
					xxx[iii], yyy[iii],
					xxx[iii+n], yyy[iii+n], dd);
 
			    /** Test corners for intersection with previous labels **/
			    label2 = labelList;
			    result = 0;
			    while ((result == 0) && (label2 != R_NilValue)) {
				result = TestLabelIntersection(label1, CAR(label2));
				label2 = CDR(label2);
			    }
			    if (result == 0) {
				result = LabelInsideWindow(label1, dd);
				if (result == 0) {
				    indx = iii;
				    range = n;
				    gotLabel = TRUE;
				}
			    }
			}
			break;
		    case 2: /* draw label on flattest portion of contour
			       embedded in contour line
			       no overlapping labels
			    */
			/* Look for flatest sequence of contour gradients */
			lowestVariance = 9999999;   /* A large number */
			indx = 0;
			range = 0;
			gotLabel = FALSE;
			for (iii = 0; iii < ns; iii++) {
			    distanceSum = 0;
			    avgGradient = 0;
			    squareSum = 0;
			    sum = 0;
			    n = 0;
			    zeroCount = 0;
			    jjj = (iii + 1);
			    while ((jjj < ns-1) &&
				   (distanceSum < labelDistance)) {
 
				/* Find a gap big enough for the label
				   use several segments if necessary
				*/
				dX = xxx[jjj] - xxx[jjj - n - 1];
				dY = yyy[jjj] - yyy[jjj - n - 1];
				dXC = GConvertXUnits(dX, USER, INCHES, dd);
				dYC = GConvertYUnits(dY, USER, INCHES, dd);
				distanceSum = hypot(dXC, dYC);
 
				/* Calculate the variance of the gradients
				   of the segments that will make way for the
				   label
				*/
				deltaX = xxx[jjj] - xxx[jjj - 1];
				deltaY = yyy[jjj] - yyy[jjj - 1];
				deltaXC = GConvertXUnits(deltaX, USER, INCHES, dd);
				deltaYC = GConvertYUnits(deltaY, USER, INCHES, dd);
				if (deltaX == 0) {deltaX = 1;}
				avgGradient += (deltaY/deltaX);
				squareSum += avgGradient * avgGradient;
				jjj = (jjj + 1);
				n += 1;
			    }
			    if (distanceSum < labelDistance)
				break;
 
			    /** Find 4 corners of label extents **/
			    FindCorners(labelDistance, labelHeight, label1,
					xxx[iii], yyy[iii],
					xxx[iii+n], yyy[iii+n], dd);
 
			    /** Test corners for intersection with previous labels **/
			    label2 = labelList;
			    result = 0;
			    while ((result == 0) && (label2 != R_NilValue)) {
				result = TestLabelIntersection(label1, CAR(label2));
				label2 = CDR(label2);
			    }
			    if (result == 0)
				result = LabelInsideWindow(label1, dd);
			    if (result == 0) {
				variance = (squareSum - (avgGradient * avgGradient) / n) / n;
				avgGradient /= n;
				if (variance < lowestVariance) {
				    lowestVariance = variance;
				    indx = iii;
				    range = n;
				    avg = avgGradient;
				}
			    }
			    if (lowestVariance < 9999999)
				gotLabel = TRUE;
			}
		    } /* switch (method) */
 
		    if (method == 0) {
			GPolyline(ns, xxx, yyy, USER, dd);
			if (vectorFonts)
			    GVText(xxx[indx], yyy[indx], USER, buffer,
				   typeface, fontindex,
				   .5, .5, 0, dd);
			else
			    GText(xxx[indx], yyy[indx], USER, buffer,
				  .5, .5, 0, dd);
		    }
		    else {
			for (iii = 0; iii < indx; iii++)
			    GLine(xxx[iii], yyy[iii],
				  xxx[iii+1], yyy[iii+1], USER, dd);
			for (iii = indx+range; iii < ns - 1; iii++)
			    GLine(xxx[iii], yyy[iii],
				  xxx[iii+1], yyy[iii+1], USER, dd);
 
			if (gotLabel) {
			    /* find which plot edge we are closest to */
			    int closest; /* 0 = indx,  1 = indx+range */
			    double dx1, dx2, dy1, dy2, dmin;
			    dx1 = fmin2((xxx[indx] - Rf_gpptr(dd)->usr[0]),
					(Rf_gpptr(dd)->usr[1] - xxx[indx]));
			    dx2 = fmin2((Rf_gpptr(dd)->usr[1] - xxx[indx+range]),
					(xxx[indx+range] - Rf_gpptr(dd)->usr[0]));
			    if (dx1 < dx2) {
				closest = 0;
				dmin = dx1;
			    } else {
				closest = 1;
				dmin = dx2;
			    }
			    dy1 = fmin2((yyy[indx] - Rf_gpptr(dd)->usr[2]),
					(Rf_gpptr(dd)->usr[3] - yyy[indx]));
			    if (closest && (dy1 < dmin)) {
				closest = 0;
				dmin = dy1;
			    } else if (dy1 < dmin)
				dmin = dy1;
			    dy2 = fmin2((Rf_gpptr(dd)->usr[3] - yyy[indx+range]),
					(yyy[indx+range] - Rf_gpptr(dd)->usr[2]));
			    if (!closest && (dy2 < dmin))
				closest = 1;
 
			    dx = GConvertXUnits(xxx[indx+range] - xxx[indx],
						USER, INCHES, dd);
			    dy = GConvertYUnits(yyy[indx+range] - yyy[indx],
						USER, INCHES, dd);
			    dxy = hypot(dx, dy);
 
			    /* save the current label for checking overlap */
			    label2 = allocVector(REALSXP, 8);
 
			    FindCorners(labelDistance, labelHeight, label2,
					xxx[indx], yyy[indx],
					xxx[indx+range], yyy[indx+range], dd);
			    UNPROTECT_PTR(labelList);
			    labelList = PROTECT(CONS(label2, labelList));
 
			    ddl = FALSE;
			    /* draw an extra bit of segment if the label
			       doesn't fill the gap */
			    if (closest) {
				xStart = xxx[indx+range] -
				    (xxx[indx+range] - xxx[indx]) *
				    labelDistance / dxy;
				yStart = yyy[indx+range] -
				    (yyy[indx+range] - yyy[indx]) *
				    labelDistance / dxy;
				if (labelDistance / dxy < 1)
				    GLine(xxx[indx], yyy[indx],
					  xStart, yStart,
					  USER, dd);
			    } else {
				xStart = xxx[indx] +
				    (xxx[indx+range] - xxx[indx]) *
				    labelDistance / dxy;
				yStart = yyy[indx] +
				    (yyy[indx+range] - yyy[indx]) *
				    labelDistance / dxy;
				if (labelDistance / dxy < 1)
				    GLine(xStart, yStart,
					  xxx[indx+range], yyy[indx+range],
					  USER, dd);
			    }
 
			    /*** Draw contour labels ***/
			    if (xxx[indx] < xxx[indx+range]) {
				if (closest) {
				    ux = xStart;
				    uy = yStart;
				    vx = xxx[indx+range];
				    vy = yyy[indx+range];
				} else {
				    ux = xxx[indx];
				    uy = yyy[indx];
				    vx = xStart;
				    vy = yStart;
				}
			    }
			    else {
				if (closest) {
				    ux = xxx[indx+range];
				    uy = yyy[indx+range];
				    vx = xStart;
				    vy = yStart;
				} else {
				    ux = xStart;
				    uy = yStart;
				    vx = xxx[indx];
				    vy = yyy[indx];
				}
			    }
 
			    if (!ddl) {
				/* convert to INCHES for calculation of
				   angle to draw text
				*/
				GConvert(&ux, &uy, USER, INCHES, dd);
				GConvert(&vx, &vy, USER, INCHES, dd);
				/* 1, 1 => sans serif, basic font
				   0, .5 => left, centre justified */
				if (vectorFonts)
				    GVText(ux, uy, INCHES, buffer,
					   typeface, fontindex, 0, .5,
					   (180 / 3.14) * atan2(vy - uy, vx - ux),
					   dd);
				else
				    GText (ux, uy, INCHES, buffer, 0, .5,
					   (180 / 3.14) * atan2(vy - uy, vx - ux),
					   dd);
			    }
			} /* if (gotLabel) */
		    } /* if (method == 0) else ... */
		} /* if (labelDistance > 0) */
 
	    } /* if (drawLabels) */
	    else {
		GPolyline(ns, xxx, yyy, USER, dd);
	    }
 
	    GMode(0, dd);
	    vmaxset(vmax);
	} /* while */
      } /* for(i .. )  for(j ..) */
    vmaxset(vmax); /* now we are done with ctr_SegDB */
    UNPROTECT_PTR(label1); /* pwwwargh! This is messy, but last thing
			      protected is likely labelList, and that needs
			      to be preserved across calls */
}
 
/* contour(x, y, z, levels, labels, labcex, drawlabels, 
 *         method, vfont, col = col, lty = lty, lwd = lwd)
 */
SEXP do_contour(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP oargs, c, x, y, z, vfont, col, rawcol, lty, lwd, labels;
    int i, j, nx, ny, nc, ncol, nlty, nlwd;
    int ltysave, colsave, lwdsave;
    double cexsave;
    double atom, zmin, zmax;
    char *vmax, *vmax0;
    int method;
    Rboolean drawLabels;
    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");
 
    oargs = args;
 
    x = CAR(args);
    internalTypeCheck(call, x, REALSXP);
    nx = LENGTH(x);
    args = CDR(args);
 
    y = CAR(args);
    internalTypeCheck(call, y, REALSXP);
    ny = LENGTH(y);
    args = CDR(args);
 
    z = CAR(args);
    internalTypeCheck(call, z, REALSXP);
    args = CDR(args);
 
    /* levels */
    c = CAR(args);
    internalTypeCheck(call, c, REALSXP);
    nc = LENGTH(c);
    args = CDR(args);
 
    labels = CAR(args);
    if (!isNull(labels))
	internalTypeCheck(call, labels, STRSXP);
    args = CDR(args);
 
    labcex = asReal(CAR(args));
    args = CDR(args);
 
    drawLabels = (Rboolean)asLogical(CAR(args));
    args = CDR(args);
 
    method = asInteger(CAR(args)); args = CDR(args);
    if (method < 1 || method > 3)
	errorcall(call, "invalid value for \"method\"");
 
    PROTECT(vfont = FixupVFont(CAR(args)));
    if (!isNull(vfont)) {
	vectorFonts = TRUE;
	typeface = INTEGER(vfont)[0];
	fontindex = INTEGER(vfont)[1];
    }
    args = CDR(args);
 
    rawcol = CAR(args);
    PROTECT(col = FixupCol(rawcol, NA_INTEGER));
    ncol = length(col);
    args = CDR(args);
 
    PROTECT(lty = FixupLty(CAR(args), Rf_gpptr(dd)->lty));
    nlty = length(lty);
    args = CDR(args);
 
    PROTECT(lwd = FixupLwd(CAR(args), Rf_gpptr(dd)->lwd));
    nlwd = length(lwd);
    args = CDR(args);
 
    if (nx < 2 || ny < 2)
	errorcall(call, "insufficient x or y values");
 
    if (nrows(z) != nx || ncols(z) != ny)
	errorcall(call, "dimension mismatch");
 
    if (nc < 1)
	errorcall(call, "no contour values");
 
    for (i = 0; i < nx; i++) {
	if (!R_FINITE(REAL(x)[i]))
	    errorcall(call, "missing x values");
	if (i > 0 && REAL(x)[i] < REAL(x)[i - 1])
	    errorcall(call, "increasing x values expected");
    }
 
    for (i = 0; i < ny; i++) {
	if (!R_FINITE(REAL(y)[i]))
	    errorcall(call, "missing y values");
	if (i > 0 && REAL(y)[i] < REAL(y)[i - 1])
	    errorcall(call, "increasing y values expected");
    }
 
    for (i = 0; i < nc; i++)
	if (!R_FINITE(REAL(c)[i]))
	    errorcall(call, "illegal NA contour values");
 
    zmin = DBL_MAX;
    zmax = DBL_MIN;
    for (i = 0; i < nx * ny; i++)
	if (R_FINITE(REAL(z)[i])) {
	    if (zmax < REAL(z)[i]) zmax =  REAL(z)[i];
	    if (zmin > REAL(z)[i]) zmin =  REAL(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);
 
    /* 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 */
 
    vmax0 = vmaxget();
    ctr_SegDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));
 
    for (i = 0; i < nx; i++)
	for (j = 0; j < ny; j++)
	    ctr_SegDB[i + j * nx] = NULL;
 
    /* Draw the contours -- note the heap release */
 
    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++) {
	vmax = vmaxget();
	Rf_gpptr(dd)->lty = INTEGER(lty)[i % nlty];
	if (Rf_gpptr(dd)->lty == NA_INTEGER)
	    Rf_gpptr(dd)->lty = ltysave;
 
	if (isNAcol(rawcol, i, ncol))
	    Rf_gpptr(dd)->col = colsave;
	else 
	    Rf_gpptr(dd)->col = INTEGER(col)[i % ncol];
	Rf_gpptr(dd)->lwd = REAL(lwd)[i % nlwd];
	if (Rf_gpptr(dd)->lwd == NA_REAL)
	    Rf_gpptr(dd)->lwd = lwdsave;
	Rf_gpptr(dd)->cex = labcex;
	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 result;
}
 
 
	/*  F i l l e d   C o n t o u r   P l o t s  */
 
	/*  R o s s  I h a k a,  M a r c h  1 9 9 9  */
 
static void
FindCutPoints(double low, double high,
	      double x1, double y1, double z1,
	      double x2, double y2, double z2,
	      double *x, double *y, double *z,
	      int *npt)
{
    double c;
 
    if (z1 > z2 ) {
	if (z2 > high || z1 < low)
	    return;
	if (z1 < high) {
	    x[*npt] = x1;
	    y[*npt] = y1;
	    z[*npt] = z1;
	    ++*npt;
	}
	else {
	    c = (z1 - high) / (z1 - z2);
	    x[*npt] = x1 + c * (x2 - x1);
	    y[*npt] = y1;
	    z[*npt] = z1 + c * (z2 - z1);
	    ++*npt;
	}
	if (z2 > low) {
	}
	else {
	    c = (z2 -low) / (z2 - z1);
	    x[*npt] = x2 - c * (x2 - x1);
	    y[*npt] = y1;
	    z[*npt] = z2 - c * (z2 - z1);
	    ++*npt;
	}
    }
    else if (z1 < z2) {
	if (z2 < low || z1 > high)
	    return;
	if (z1 > low) {
	    x[*npt] = x1;
	    y[*npt] = y1;
	    z[*npt] = z1;
	    ++*npt;
	}
	else {
	    c = (z1 - low) / (z1 - z2);
	    x[*npt] = x1 + c * (x2 - x1);
	    y[*npt] = y1;
	    z[*npt] = z1 + c * (z2 - z1);
	    ++*npt;
	}
	if (z2 < high) {
#ifdef OMIT
	    /* Don't repeat corner vertices */
	    x[*npt] = x2;
	    y[*npt] = y2;
	    z[*npt] = z2;
	    ++*npt;
#endif
	}
	else {
	    c = (z2 - high) / (z2 - z1);
	    x[*npt] = x2 - c * (x2 - x1);
	    y[*npt] = y1;
	    z[*npt] = z2 - c * (z2 - z1);
	    ++*npt;
	}
    }
    else {
	if(low <= z1 && z1 <= high) {
	    x[*npt] = x1;
	    y[*npt] = y1;
	    z[*npt] = z1;
	    ++*npt;
#ifdef OMIT
	    /* Don't repeat corner vertices */
	    x[*npt] = x2;
	    y[*npt] = y2;
	    z[*npt] = z2;
	    ++*npt;
#endif
	}
    }
}
 
/* FIXME - This could pretty easily be adapted to handle NA */
/* values on the grid.  Just search the diagonals for cutpoints */
/* instead of the cell sides.  Use the same switch idea as in */
/* contour above.  There are 5 cases to handle. */
 
static void
FindPolygonVertices(double low, double high,
		    double x1, double x2, double y1, double y2,
		    double z11, double z21, double z12, double z22,
		    double *x, double *y, double *z, int *npt)
{
    *npt = 0;
    FindCutPoints(low, high, x1,  y1,  z11, x2,  y1,  z21, x, y, z, npt);
    FindCutPoints(low, high, y1,  x2,  z21, y2,  x2,  z22, y, x, z, npt);
    FindCutPoints(low, high, x2,  y2,  z22, x1,  y2,  z12, x, y, z, npt);
    FindCutPoints(low, high, y2,  x1,  z12, y1,  x1,  z11, y, x, z, npt);
}
 
/* FIXME: [Code consistency] Use macro for the parallel parts of
	  do_contour, do_filledcontour & do_image ...
*/
 
/* filledcontour(x, y, z, levels, col) */
SEXP do_filledcontour(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP oargs, sx, sy, sz, sc, scol;
    double *x, *y, *z, *c;
    unsigned *col;
    int i, j, k, npt, nx, ny, nz, nc, ncol, colsave, xpdsave;
    double px[8], py[8], pz[8];
    DevDesc *dd = CurrentDevice();
 
    GCheckState(dd);
 
    checkArity(op,args);
    oargs = args;
 
    sx = CAR(args);
    internalTypeCheck(call, sx, REALSXP);
    nx = LENGTH(sx);
    args = CDR(args);
 
    sy = CAR(args);
    internalTypeCheck(call, sy, REALSXP);
    ny = LENGTH(sy);
    args = CDR(args);
 
    sz = CAR(args);
    internalTypeCheck(call, sz, REALSXP);
    nz = length(sz);
    args = CDR(args);
 
    sc = CAR(args);/* levels */
    internalTypeCheck(call, sc, REALSXP);
    nc = length(sc);
    args = CDR(args);
 
    if (nx < 2 || ny < 2)
	errorcall(call, "insufficient x or y values");
 
    if (nrows(sz) != nx || ncols(sz) != ny)
	errorcall(call, "dimension mismatch");
 
    if (nc < 1)
	errorcall(call, "no contour values");
 
    PROTECT(scol = FixupCol(CAR(args), NA_INTEGER));
    ncol = length(scol);
 
    /* Shorthand Pointers */
 
    x = REAL(sx);
    y = REAL(sy);
    z = REAL(sz);
    c = REAL(sc);
    col = (unsigned*)INTEGER(scol);
 
    /* Check of grid coordinates */
    /* We want them to all be finite */
    /* and in strictly ascending order */
 
    if (nx < 1 || ny < 1) goto badxy;
    if (!R_FINITE(x[0])) goto badxy;
    if (!R_FINITE(y[0])) goto badxy;
    for (i = 1; i < nx; i++)
	if (!R_FINITE(x[i]) || x[i] <= x[i - 1]) goto badxy;
    for (j = 1; j < ny; j++)
	if (!R_FINITE(y[j]) || y[j] <= y[j - 1]) goto badxy;
 
    /* Check of the contour levels */
 
    if (!R_FINITE(c[0])) goto badlev;
    for (k = 1; k < nc; k++)
	if (!R_FINITE(c[k]) || c[k] <= c[k - 1]) goto badlev;
 
    colsave = Rf_gpptr(dd)->col;
    xpdsave = Rf_gpptr(dd)->xpd;
    /* override par("xpd") and force clipping to plot region */
    Rf_gpptr(dd)->xpd = 0;
 
    GMode(1, dd);
 
    for (i = 1; i < nx; i++) {
	for (j = 1; j < ny; j++) {
	    for (k = 1; k < nc ; k++) {
		FindPolygonVertices(c[k - 1], c[k],
				    x[i - 1], x[i],
				    y[j - 1], y[j],
				    z[i - 1 + (j - 1) * nx],
				    z[i + (j - 1) * nx],
				    z[i - 1 + j * nx],
				    z[i + j * nx],
				    px, py, pz, &npt);
                if (npt > 2)
		    GPolygon(npt, px, py, USER, col[(k-1)%ncol],
			     NA_INTEGER, dd);
	    }
	}
    }
    GMode(0, dd);
    Rf_gpptr(dd)->col = colsave;
    Rf_gpptr(dd)->xpd = xpdsave;
    R_Visible = 0;
    UNPROTECT(1);
    if (GRecording(call))
	recordGraphicOperation(op, oargs, dd);
    return R_NilValue;
 
 badxy:
    errorcall(call, "invalid x / y values or limits");
 badlev:
    errorcall(call, "invalid contour levels: must be strictly increasing");
    return R_NilValue;  /* never used; to keep -Wall happy */
}
 
 
	/*  I m a g e   R e n d e r i n g  */
 
 
/* image(x, y, z, col, breaks) */
SEXP do_image(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP oargs, sx, sy, sz, sc;
    double *x, *y;
    int *z, tmp;
    unsigned *c;
    int i, j, nx, ny, nc, colsave, xpdsave;
    DevDesc *dd = CurrentDevice();
 
    GCheckState(dd);
 
    checkArity(op,args);
    oargs = args;
 
    sx = CAR(args);
    internalTypeCheck(call, sx, REALSXP);
    nx = LENGTH(sx);
    args = CDR(args);
 
    sy = CAR(args);
    internalTypeCheck(call, sy, REALSXP);
    ny = LENGTH(sy);
    args = CDR(args);
 
    sz = CAR(args);
    internalTypeCheck(call, sz, INTSXP);
    args = CDR(args);
 
    PROTECT(sc = FixupCol(CAR(args), NA_INTEGER));
    nc = LENGTH(sc);
 
    /* Shorthand Pointers */
 
    x = REAL(sx);
    y = REAL(sy);
    z = INTEGER(sz);
    c = (unsigned*)INTEGER(sc);
 
    /* Check of grid coordinates */
    /* We want them to all be finite and in strictly ascending order */
 
    if (nx < 1 || ny < 1) goto badxy;
    if (!R_FINITE(x[0])) goto badxy;
    if (!R_FINITE(y[0])) goto badxy;
    for (i = 1; i < nx; i++)
	if (!R_FINITE(x[i]) || x[i] <= x[i - 1]) goto badxy;
    for (j = 1; j < ny; j++)
	if (!R_FINITE(y[j]) || y[j] <= y[j - 1]) goto badxy;
 
    colsave = Rf_gpptr(dd)->col;
    xpdsave = Rf_gpptr(dd)->xpd;
    /* override par("xpd") and force clipping to plot region */
    Rf_gpptr(dd)->xpd = 0;
 
    GMode(1, dd);
 
    for (i = 0; i < nx - 1 ; i++) {
	for (j = 0; j < ny - 1; j++) {
	    tmp = z[i + j * (nx - 1)];
	    if (tmp >= 0 && tmp < nc && tmp != NA_INTEGER)
		GRect(x[i], y[j], x[i+1], y[j+1], USER, c[tmp],
		      NA_INTEGER, dd);
	}
    }
    GMode(0, dd);
    Rf_gpptr(dd)->col = colsave;
    Rf_gpptr(dd)->xpd = xpdsave;
    R_Visible = 0;
    UNPROTECT(1);
    if (GRecording(call))
	recordGraphicOperation(op, oargs, dd);
    return R_NilValue;
 
  badxy:
    errorcall(call, "invalid x / y values or limits");
    return R_NilValue;/* never used; to keep -Wall happy */
}
 
	/*  P e r s p e c t i v e   S u r f a c e   P l o t s  */
 
 
/* Set up the light source */
static double Light[4];
static double Shade;
static Rboolean DoLighting;
 
static void SetUpLight(double theta, double phi)
{
    double u[4];
    u[0] = 0; u[1] = -1; u[2] = 0; u[3] = 1;
    SetToIdentity(VT);             /* Initialization */
    XRotate(-phi);                 /* colatitude rotation */
    ZRotate(theta);                /* azimuthal rotation */
    TransVector(u, VT, Light);	   /* transform */
}
 
static double FacetShade(double *u, double *v)
{
    double nx, ny, nz, sum;
    nx = u[1] * v[2] - u[2] * v[1];
    ny = u[2] * v[0] - u[0] * v[2];
    nz = u[0] * v[1] - u[1] * v[0];
    sum = sqrt(nx * nx + ny * ny + nz * nz);
    if (sum == 0) sum = 1;
    nx /= sum;
    ny /= sum;
    nz /= sum;
    sum = 0.5 * (nx * Light[0] + ny * Light[1] + nz * Light[2] + 1);
    return pow(sum, Shade);
}
 
 
/* For each facet, determine the farthest point from the eye. */
/* Sorting the facets so that these depths are decreasing */
/* yields an occlusion compatible ordering. */
/* Note that we ignore z values when doing this. */
 
static void DepthOrder(double *z, double *x, double *y, int nx, int ny,
		       double *depth, int *indx)
{
    int i, ii, j, jj, nx1, ny1;
    Vector3d u, v;
    double d;
    nx1 = nx - 1;
    ny1 = ny - 1;
    for (i = 0; i < nx1 * ny1; i++)
	indx[i] = i;
    for (i = 0; i < nx1; i++)
	for (j = 0; j < ny1; j++) {
	    d = -DBL_MAX;
	    for (ii = 0; ii <= 1; ii++)
		for (jj = 0; jj <= 1; jj++) {
		    u[0] = x[i + ii];
		    u[1] = y[j + jj];
		    /* Originally I had the following line here: */
		    /* u[2] = z[i+ii+(j+jj)*nx]; */
		    /* But this leads to artifacts. */
		    /* It has been replaced by the following line: */
		    u[2] = 0;
		    u[3] = 1;
		    if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
			TransVector(u, VT, v);
			if (v[3] > d) d = v[3];
		    }
		}
	    depth[i+j*nx1] = -d;
 
	}
    /* Determine the depth ordering of the facets to ensure
       that they are drawn in an occlusion compatible order. */
    rsort_with_index(depth, indx, nx1 * ny1);
}
 
 
static void DrawFacets(double *z, double *x, double *y, int nx, int ny,
		       int *indx, double xs, double ys, double zs,
	               int *col, int ncol, int border)
{
    double xx[4], yy[4], shade = 0;
    Vector3d u, v;
    int i, j, k, n, nx1, ny1, icol, nv;
    unsigned int newcol, r, g, b;
    DevDesc *dd;
    dd = CurrentDevice();
    nx1 = nx - 1;
    ny1 = ny - 1;
    n = nx1 * ny1;
    for (k = 0; k < n; k++) {
	nv = 0;
	i = indx[k] % nx1;
	j = indx[k] / nx1;
	icol = (i + j * nx1) % ncol;
	if (DoLighting) {
            /* Note we must scale here */
	    u[0] = xs * (x[i+1] - x[i]);
	    u[1] = ys * (y[j] - y[j+1]);
	    u[2] = zs * (z[(i+1)+j*nx] - z[i+(j+1)*nx]);
	    v[0] = xs * (x[i+1] - x[i]);
	    v[1] = ys * (y[j+1] - y[j]);
	    v[2] = zs * (z[(i+1)+(j+1)*nx] - z[i+j*nx]);
	    shade = FacetShade(u, v);
	}
	u[0] = x[i]; u[1] = y[j];
	u[2] = z[i + j * nx]; u[3] = 1;
	if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
	    TransVector(u, VT, v);
	    xx[nv] = v[0] / v[3];
	    yy[nv] = v[1] / v[3];
	    nv++;
	}
 
	u[0] = x[i + 1]; u[1] = y[j];
	u[2] = z[i + 1 + j * nx]; u[3] = 1;
	if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
	    TransVector(u, VT, v);
	    xx[nv] = v[0] / v[3];
	    yy[nv] = v[1] / v[3];
	    nv++;
	}
 
	u[0] = x[i + 1]; u[1] = y[j + 1];
	u[2] = z[i + 1 + (j + 1) * nx]; u[3] = 1;
	if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
	    TransVector(u, VT, v);
	    xx[nv] = v[0] / v[3];
	    yy[nv] = v[1] / v[3];
	    nv++;
	}
 
	u[0] = x[i]; u[1] = y[j + 1];
	u[2] = z[i + (j + 1) * nx]; u[3] = 1;
	if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
	    TransVector(u, VT, v);
	    xx[nv] = v[0] / v[3];
	    yy[nv] = v[1] / v[3];
	    nv++;
	}
 
	if (nv > 2) {
	    newcol = col[icol];
	    if (DoLighting) {
		r = shade * R_RED(newcol);
		g = shade * R_GREEN(newcol);
		b = shade * R_BLUE(newcol);
		newcol = R_RGB(r, g, b);
	    }
	    GPolygon(nv, xx, yy, USER, newcol, border, dd);
	}
    }
}
 
 
#ifdef NOT_used_currently/*-- out 'def'  (-Wall) --*/
static void CheckRange(double *x, int n, double min, double max)
{
    double xmin, xmax;
    int i;
    xmin =  DBL_MAX;
    xmax = -DBL_MAX;
    for (i = 0; i < n; i++)
	if (R_FINITE(x[i])) {
	    if (x[i] < xmin) xmin = x[i];
	    if (x[i] > xmax) xmax = x[i];
	}
    if (xmin < min || xmax > max)
	errorcall(gcall, "coordinates outsize specified range");
}
#endif
 
static void PerspWindow(double *xlim, double *ylim, double *zlim, DevDesc *dd)
{
    double pin1, pin2, scale, xdelta, ydelta, xscale, yscale, xadd, yadd;
    double xmax, xmin, ymax, ymin, xx, yy;
    Vector3d u, v;
    int i, j, k;
 
    xmax = xmin = ymax = ymin = 0;
    u[3] = 1;
    for (i = 0; i < 2; i++) {
	u[0] = xlim[i];
	for (j = 0; j < 2; j++) {
	    u[1] = ylim[j];
	    for (k = 0; k < 2; k++) {
		u[2] = zlim[k];
		TransVector(u, VT, v);
		xx = v[0] / v[3];
		yy = v[1] / v[3];
		if (xx > xmax) xmax = xx;
		if (xx < xmin) xmin = xx;
		if (yy > ymax) ymax = yy;
		if (yy < ymin) ymin = yy;
	    }
	}
    }
    pin1 = GConvertXUnits(1.0, NPC, INCHES, dd);
    pin2 = GConvertYUnits(1.0, NPC, INCHES, dd);
    xdelta = fabs(xmax - xmin);
    ydelta = fabs(ymax - ymin);
    xscale = pin1 / xdelta;
    yscale = pin2 / ydelta;
    scale = (xscale < yscale) ? xscale : yscale;
    xadd = .5 * (pin1 / scale - xdelta);
    yadd = .5 * (pin2 / scale - ydelta);
    GScale(xmin - xadd, xmax + xadd, 1, dd);
    GScale(ymin - yadd, ymax + yadd, 2, dd);
    GMapWin2Fig(dd);
}
 
static int LimitCheck(double *lim, double *c, double *s)
{
    if (!R_FINITE(lim[0]) || !R_FINITE(lim[1]) || lim[0] >= lim[1])
	return 0;
    *s = 0.5 * fabs(lim[1] - lim[0]);
    *c = 0.5 * (lim[1] + lim[0]);
    return 1;
}
 
/* PerspBox: The following code carries out a visibility test */
/* on the surfaces of the xlim/ylim/zlim box around the plot. */
/* If front = 0, only the faces with their inside toward the */
/* eyepoint are drawn.  If front = 1, only the faces with */
/* their outside toward the eye are drawn.  This lets us carry */
/* out hidden line removal by drawing any faces which will be */
/* obscured before the surface, and those which will not be */
/* obscured after the surface. */
 
/* The vertices of the box */
static short int Vertex[8][3] = {
    {0, 0, 0},
    {0, 0, 1},
    {0, 1, 0},
    {0, 1, 1},
    {1, 0, 0},
    {1, 0, 1},
    {1, 1, 0},
    {1, 1, 1},
};
 
/* The vertices visited when tracing a face */
static short int Face[6][4] = {
    {0, 1, 5, 4},
    {2, 6, 7, 3},
    {0, 2, 3, 1},
    {4, 5, 7, 6},
    {0, 4, 6, 2},
    {1, 3, 7, 5},
};
 
/* The edges drawn when tracing a face */
static short int Edge[6][4] = {
    { 0, 1, 2, 3},
    { 4, 5, 6, 7},
    { 8, 7, 9, 0},
    { 2,10, 5,11},
    { 3,11, 4, 8},
    { 9, 6,10, 1},
};
 
/* Which edges have been drawn previously */
static char EdgeDone[12];
 
static void PerspBox(int front, double *x, double *y, double *z, DevDesc *dd)
{
    Vector3d u0, v0, u1, v1, u2, v2, u3, v3;
    double d[3], e[3];
    int f, i, p0, p1, p2, p3, nearby;
    int ltysave = Rf_gpptr(dd)->lty;
    if (front)
	Rf_gpptr(dd)->lty = LTY_DOTTED;
    else
	Rf_gpptr(dd)->lty = LTY_SOLID;
    for (f = 0; f < 6; f++) {
        p0 = Face[f][0];
        p1 = Face[f][1];
        p2 = Face[f][2];
        p3 = Face[f][3];
 
	u0[0] = x[Vertex[p0][0]];
	u0[1] = y[Vertex[p0][1]];
	u0[2] = z[Vertex[p0][2]];
	u0[3] = 1;
	u1[0] = x[Vertex[p1][0]];
	u1[1] = y[Vertex[p1][1]];
	u1[2] = z[Vertex[p1][2]];
	u1[3] = 1;
	u2[0] = x[Vertex[p2][0]];
	u2[1] = y[Vertex[p2][1]];
	u2[2] = z[Vertex[p2][2]];
	u2[3] = 1;
	u3[0] = x[Vertex[p3][0]];
	u3[1] = y[Vertex[p3][1]];
	u3[2] = z[Vertex[p3][2]];
	u3[3] = 1;
 
	TransVector(u0, VT, v0);
	TransVector(u1, VT, v1);
	TransVector(u2, VT, v2);
	TransVector(u3, VT, v3);
 
	/* Visibility test. */
	/* Determine whether the surface normal is toward the eye. */
	/* Note that we only draw lines once. */
 
        for (i = 0; i < 3; i++) {
	    d[i] = v1[i]/v1[3] - v0[i]/v0[3];
	    e[i] = v2[i]/v2[3] - v1[i]/v1[3];
        }
	nearby = (d[0]*e[1] - d[1]*e[0]) < 0;
 
	if ((front && nearby) || (!front && !nearby)) {
	    if (!EdgeDone[Edge[f][0]]++)
		GLine(v0[0]/v0[3], v0[1]/v0[3],
		      v1[0]/v1[3], v1[1]/v1[3], USER, dd);
	    if (!EdgeDone[Edge[f][1]]++)
		GLine(v1[0]/v1[3], v1[1]/v1[3],
		      v2[0]/v2[3], v2[1]/v2[3], USER, dd);
	    if (!EdgeDone[Edge[f][2]]++)
		GLine(v2[0]/v2[3], v2[1]/v2[3],
		      v3[0]/v3[3], v3[1]/v3[3], USER, dd);
	    if (!EdgeDone[Edge[f][3]]++)
		GLine(v3[0]/v3[3], v3[1]/v3[3],
		      v0[0]/v0[3], v0[1]/v0[3], USER, dd);
	}
    }
    Rf_gpptr(dd)->lty = ltysave;
}
 
/* PerspAxes:
 */
 
/* Starting vertex for possible axes */
static short int AxisStart[8] = { 0, 0, 2, 4, 0, 4, 2, 6 };
 
/* Tick vector for possible axes */
static short int TickVector[8][3] = {
    {0, -1, -1},
    {-1, 0, -1},
    {0, 1, -1},
    {1, 0, -1},
    {-1, -1, 0},
    {1, -1, 0},
    {-1, 1, 0},
    {1, 1, 0}};
 
static int lowest(double y1, double y2, double y3, double y4) {
    return ((y1 <= y2) && (y1 <= y3) && (y1 <= y4));
}
 
static double labelAngle(double x1, double y1, double x2, double y2) {
    double dx, dy;
    double angle;
    dx = fabs(x2 - x1);
    if (x2 > x1)
	dy = y2 - y1;
    else
	dy = y1 - y2;
    if (dx == 0) {
	if (dy > 0)
	    angle = 90;
	else
	    angle = 270;
    } else {
	angle = (180 / M_PI) * atan2(dy, dx);
    }
    return angle;
}
 
static void PerspAxis(double *x, double *y, double *z,
		      int axis, int axisType, int nTicks, int tickType,
		      char *label, DevDesc *dd) {
    Vector3d u1, u2, u3, v1, v2, v3;
    double tickLength = .03; /* proportion of axis length */
    double min, max, d_frac;
    double *range = NULL; /* -Wall */
    double axp[3];
    int nint, i;
    SEXP at, lab;
    switch (axisType) {
    case 0:
	min = x[0];	max = x[1];	range = x;	break;
    case 1:
	min = y[0];	max = y[1];	range = y;	break;
    case 2:
	min = z[0];	max = z[1];	range = z;	break;
    }
    d_frac = 0.1*(max - min);
    nint = nTicks - 1; if(!nint) nint++;
    i = nint;
    GPretty(&min, &max, &nint);
    /* GPretty() rarely gives values too much outside range ..
       2D axis() clip these, we play cheaper */ 
    while((min < range[0] - d_frac || range[1] + d_frac < max) && i < 20) {
	nint = ++i;
	min = range[0];
	max = range[1];
	GPretty(&min, &max, &nint);
    }
    axp[0] = min;
    axp[1] = max;
    axp[2] = nint;
    /* Do the following calculations for both ticktypes */
    switch (axisType) {
    case 0:
	u1[0] = min;
	u1[1] = y[Vertex[AxisStart[axis]][1]];
	u1[2] = z[Vertex[AxisStart[axis]][2]];
	break;
    case 1:
	u1[0] = x[Vertex[AxisStart[axis]][0]];
	u1[1] = min;
	u1[2] = z[Vertex[AxisStart[axis]][2]];
	break;
    case 2:
	u1[0] = x[Vertex[AxisStart[axis]][0]];
	u1[1] = y[Vertex[AxisStart[axis]][1]];
	u1[2] = min;
	break;
    }
    u1[0] = u1[0] + tickLength*(x[1]-x[0])*TickVector[axis][0];
    u1[1] = u1[1] + tickLength*(y[1]-y[0])*TickVector[axis][1];
    u1[2] = u1[2] + tickLength*(z[1]-z[0])*TickVector[axis][2];
    u1[3] = 1;
    switch (axisType) {
    case 0:
	u2[0] = max;
	u2[1] = u1[1];
	u2[2] = u1[2];
	break;
    case 1:
	u2[0] = u1[0];
	u2[1] = max;
	u2[2] = u1[2];
	break;
    case 2:
	u2[0] = u1[0];
	u2[1] = u1[1];
	u2[2] = max;
	break;
    }
    u2[3] = 1;
    /* The axis label has to be further out for "detailed" ticks
       in order to leave room for the tick labels */
    switch (tickType) {
    case 1: /* "simple": just an arrow parallel to axis, indicating direction
	       of increase */
	u3[0] = u1[0] + tickLength*(x[1]-x[0])*TickVector[axis][0];
	u3[1] = u1[1] + tickLength*(y[1]-y[0])*TickVector[axis][1];
	u3[2] = u1[2] + tickLength*(z[1]-z[0])*TickVector[axis][2];
	break;
    case 2:
	u3[0] = u1[0] + 2.5*tickLength*(x[1]-x[0])*TickVector[axis][0];
	u3[1] = u1[1] + 2.5*tickLength*(y[1]-y[0])*TickVector[axis][1];
	u3[2] = u1[2] + 2.5*tickLength*(z[1]-z[0])*TickVector[axis][2];
	break;
    }
    switch (axisType) {
    case 0:
	u3[0] = (min + max)/2;
	break;
    case 1:
	u3[1] = (min + max)/2;
	break;
    case 2:
	u3[2] = (min + max)/2;
	break;
    }
    u3[3] = 1;
    TransVector(u1, VT, v1);
    TransVector(u2, VT, v2);
    TransVector(u3, VT, v3);
    /* Draw axis label */
    GText(v3[0]/v3[3], v3[1]/v3[3], USER, label, .5, .5,
	  labelAngle(v1[0]/v1[3], v1[1]/v1[3], v2[0]/v2[3], v2[1]/v2[3]),
	  dd);
    /* Draw axis ticks */
    switch (tickType) {
    case 1: /* "simple": just an arrow parallel to axis, indicating direction
	       of increase */
	/* arrow head is 0.25 inches long, with angle 30 degrees,
	   and drawn at v2 end of line */
	GArrow(v1[0]/v1[3], v1[1]/v1[3],
	       v2[0]/v2[3], v2[1]/v2[3], USER,
	       0.1, 10, 2, dd);
	break;
    case 2: /* "detailed": normal ticks as per 2D plots */
	PROTECT(at = CreateAtVector(axp, range, 7, FALSE));
	PROTECT(lab = labelformat(at));
	for (i=0; i<length(at); i++) {
	    switch (axisType) {
	    case 0:
		u1[0] = REAL(at)[i];
		u1[1] = y[Vertex[AxisStart[axis]][1]];
		u1[2] = z[Vertex[AxisStart[axis]][2]];
		break;
	    case 1:
		u1[0] = x[Vertex[AxisStart[axis]][0]];
		u1[1] = REAL(at)[i];
		u1[2] = z[Vertex[AxisStart[axis]][2]];
		break;
	    case 2:
		u1[0] = x[Vertex[AxisStart[axis]][0]];
		u1[1] = y[Vertex[AxisStart[axis]][1]];
		u1[2] = REAL(at)[i];
		break;
	    }
	    u1[3] = 1;
	    u2[0] = u1[0] + tickLength*(x[1]-x[0])*TickVector[axis][0];
	    u2[1] = u1[1] + tickLength*(y[1]-y[0])*TickVector[axis][1];
	    u2[2] = u1[2] + tickLength*(z[1]-z[0])*TickVector[axis][2];
	    u2[3] = 1;
	    u3[0] = u2[0] + tickLength*(x[1]-x[0])*TickVector[axis][0];
	    u3[1] = u2[1] + tickLength*(y[1]-y[0])*TickVector[axis][1];
	    u3[2] = u2[2] + tickLength*(z[1]-z[0])*TickVector[axis][2];
	    u3[3] = 1;
	    TransVector(u1, VT, v1);
	    TransVector(u2, VT, v2);
	    TransVector(u3, VT, v3);
	    /* Draw tick line */
	    GLine(v1[0]/v1[3], v1[1]/v1[3],
		  v2[0]/v2[3], v2[1]/v2[3], USER, dd);
	    /* Draw tick label */
	    GText(v3[0]/v3[3], v3[1]/v3[3], USER, CHAR(STRING_ELT(lab, i)),
		  .5, .5, 0, dd);
	}
	UNPROTECT(2);
	break;
    }
}
 
/* Determine the transformed (x, y) coordinates (in USER space)
 * for the four corners of the x-y plane of the persp plot
 * These will be used to determine which sides of the persp
 * plot to label with axes
 * The strategy is to determine which corner has the lowest y-value
 * to decide which of the x- and y-axes to label AND which corner
 * has the lowest x-value to decide which of the z-axes to label
 */
static void PerspAxes(double *x, double *y, double *z,
                      char *xlab, char *ylab, char *zlab,
		      int nTicks, int tickType, DevDesc *dd) {
    int xAxis=0, yAxis=0, zAxis=0; /* -Wall */
    int xpdsave;
    Vector3d u0, u1, u2, u3;
    Vector3d v0, v1, v2, v3;
    u0[0] = x[0];
    u0[1] = y[0];
    u0[2] = z[0];
    u0[3] = 1;
    u1[0] = x[1];
    u1[1] = y[0];
    u1[2] = z[0];
    u1[3] = 1;
    u2[0] = x[0];
    u2[1] = y[1];
    u2[2] = z[0];
    u2[3] = 1;
    u3[0] = x[1];
    u3[1] = y[1];
    u3[2] = z[0];
    u3[3] = 1;
    TransVector(u0, VT, v0);
    TransVector(u1, VT, v1);
    TransVector(u2, VT, v2);
    TransVector(u3, VT, v3);
 
    /* to fit in the axis labels */
    xpdsave = Rf_gpptr(dd)->xpd;
    Rf_gpptr(dd)->xpd = 1;
 
    /* Figure out which X and Y axis to draw */
    if (lowest(v0[1]/v0[3], v1[1]/v1[3], v2[1]/v2[3], v3[1]/v3[3])) {
	xAxis = 0;
	yAxis = 1;
    } else if (lowest(v1[1]/v1[3], v0[1]/v0[3], v2[1]/v2[3], v3[1]/v3[3])) {
	xAxis = 0;
	yAxis = 3;
    } else if (lowest(v2[1]/v2[3], v1[1]/v1[3], v0[1]/v0[3], v3[1]/v3[3])) {
	xAxis = 2;
	yAxis = 1;
    } else if (lowest(v3[1]/v3[3], v1[1]/v1[3], v2[1]/v2[3], v0[1]/v0[3])) {
	xAxis = 2;
	yAxis = 3;
    } else
	warning("Axis orientation not calculated");
    PerspAxis(x, y, z, xAxis, 0, nTicks, tickType, xlab, dd);
    PerspAxis(x, y, z, yAxis, 1, nTicks, tickType, ylab, dd);
    /* Figure out which Z axis to draw */
    if (lowest(v0[0]/v0[3], v1[0]/v1[3], v2[0]/v2[3], v3[0]/v3[3])) {
	zAxis = 4;
    } else if (lowest(v1[0]/v1[3], v0[0]/v0[3], v2[0]/v2[3], v3[0]/v3[3])) {
	zAxis = 5;
    } else if (lowest(v2[0]/v2[3], v1[0]/v1[3], v0[0]/v0[3], v3[0]/v3[3])) {
	zAxis = 6;
    } else if (lowest(v3[0]/v3[3], v1[0]/v1[3], v2[0]/v2[3], v0[0]/v0[3])) {
	zAxis = 7;
    } else
	warning("Axes orientation not calculated");
    PerspAxis(x, y, z, zAxis, 2, nTicks, tickType, zlab, dd);
 
    Rf_gpptr(dd)->xpd = xpdsave;
}
 
SEXP do_persp(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP x, y, z, xlim, ylim, zlim;
    SEXP depth, indx, originalArgs;
    SEXP col, border, xlab, ylab, zlab;
    double theta, phi, r, d;
    double ltheta, lphi;
    double expand, xc, yc, zc, xs, ys, zs;
    int i, j, scale, ncol, dobox, doaxes, nTicks, tickType;
    DevDesc *dd;
 
    if (length(args) < 24)
	errorcall(call, "too few parameters");
    gcall = call;
    originalArgs = args;
 
    PROTECT(x = coerceVector(CAR(args), REALSXP));
    if (length(x) < 2) errorcall(call, "invalid x argument");
    args = CDR(args);
 
    PROTECT(y = coerceVector(CAR(args), REALSXP));
    if (length(y) < 2) errorcall(call, "invalid y argument");
    args = CDR(args);
 
    PROTECT(z = coerceVector(CAR(args), REALSXP));
    if (!isMatrix(z) || nrows(z) != length(x) || ncols(z) != length(y))
	errorcall(call, "invalid z argument");
    args = CDR(args);
 
    PROTECT(xlim = coerceVector(CAR(args), REALSXP));
    if (length(xlim) != 2) errorcall(call, "invalid xlim argument");
    args = CDR(args);
 
    PROTECT(ylim = coerceVector(CAR(args), REALSXP));
    if (length(ylim) != 2) errorcall(call, "invalid ylim argument");
    args = CDR(args);
 
    PROTECT(zlim = coerceVector(CAR(args), REALSXP));
    if (length(zlim) != 2) errorcall(call, "invalid zlim argument");
    args = CDR(args);
 
    /* Checks on x/y/z Limits */
 
    if (!LimitCheck(REAL(xlim), &xc, &xs))
	errorcall(call, "invalid x limits");
    if (!LimitCheck(REAL(ylim), &yc, &ys))
	errorcall(call, "invalid y limits");
    if (!LimitCheck(REAL(zlim), &zc, &zs))
	errorcall(call, "invalid z limits");
 
    theta = asReal(CAR(args));	args = CDR(args);
    phi	  = asReal(CAR(args));	args = CDR(args);
    r	= asReal(CAR(args));	args = CDR(args);
    d	= asReal(CAR(args));	args = CDR(args);
    scale  = asLogical(CAR(args)); args = CDR(args);
    expand = asReal(CAR(args)); args = CDR(args);
    col	   = CAR(args);		args = CDR(args);
    border = CAR(args);		args = CDR(args);
    ltheta = asReal(CAR(args)); args = CDR(args);
    lphi   = asReal(CAR(args)); args = CDR(args);
    Shade  = asReal(CAR(args)); args = CDR(args);
    dobox  = asLogical(CAR(args)); args = CDR(args);
    doaxes = asLogical(CAR(args)); args = CDR(args);
    nTicks = asInteger(CAR(args)); args = CDR(args);
    tickType = asInteger(CAR(args)); args = CDR(args);
    xlab = CAR(args); args = CDR(args);
    ylab = CAR(args); args = CDR(args);
    zlab = CAR(args); args = CDR(args);
    if (!isString(xlab) || length(xlab) < 1)
	error("`xlab' must be a character vector of length 1");
    if (!isString(ylab) || length(ylab) < 1)
	error("`ylab' must be a character vector of length 1");
    if (!isString(zlab) || length(zlab) < 1)
	error("`zlab' must be a character vector of length 1");
 
    if (R_FINITE(Shade) && Shade <= 0) Shade = 1;
    if (R_FINITE(ltheta) && R_FINITE(lphi) && R_FINITE(Shade))
	DoLighting = TRUE;
    else
	DoLighting = FALSE;
 
    if (!scale) {
	double s;
	s = xs;
	if (s < ys) s = ys;
	if (s < zs) s = zs;
	xs = s; ys = s; zs = s;
    }
 
    /* Parameter Checks */
 
    if (!R_FINITE(theta) || !R_FINITE(phi) || !R_FINITE(r) || !R_FINITE(d) ||
	d < 0 || r < 0)
	errorcall(call, "invalid viewing parameters");
    if (!R_FINITE(expand) || expand < 0)
	errorcall(call, "invalid expand value");
    if (scale == NA_LOGICAL)
	scale = 0;
    if ((nTicks == NA_INTEGER) || (nTicks < 0))
	errorcall(call, "invalid nticks value");
    if ((tickType == NA_INTEGER) || (tickType < 1) || (tickType > 2))
	errorcall(call, "invalid ticktype value");
 
    dd = GNewPlot(GRecording(call));
 
    PROTECT(col = FixupCol(col, Rf_gpptr(dd)->bg));
    ncol = LENGTH(col);
    if (ncol < 1) errorcall(call, "invalid col specification");
    if(!R_OPAQUE(INTEGER(col)[0])) DoLighting = FALSE;
    PROTECT(border = FixupCol(border, Rf_gpptr(dd)->fg));
    if (length(border) < 1) errorcall(call, "invalid border specification");
 
    GSetState(1, dd);
    GSavePars(dd);
    ProcessInlinePars(args, dd, call);
    if (length(border) > 1)
	Rf_gpptr(dd)->fg = INTEGER(border)[0];
    Rf_gpptr(dd)->xlog = Rf_gpptr(dd)->ylog = FALSE;
 
    /* Set up the light vector (if any) */
    if (DoLighting)
        SetUpLight(ltheta, lphi);
 
    /* Mark box edges as undrawn */
    for (i = 0; i< 12; i++)
	EdgeDone[i] = 0;
 
    /* Specify the viewing transformation. */
 
    SetToIdentity(VT);             /* Initialization */
    Translate(-xc, -yc, -zc);      /* center at the origin */
    Scale(1/xs, 1/ys, expand/zs);  /* scale extents to [-1,1] */
    XRotate(-90.0);                /* rotate x-y plane to horizontal */
    YRotate(-theta);               /* azimuthal rotation */
    XRotate(phi);                  /* elevation rotation */
    Translate(0.0, 0.0, -r - d);   /* translate the eyepoint to the origin */
    Perspective(d);                /* perspective */
 
    /* Specify the plotting window. */
    /* Here we map the vertices of the cube */
    /* [xmin,xmax]*[ymin,ymax]*[zmin,zmax] */
    /* to the screen and then chose a window */
    /* which is symmetric about (0,0). */
 
    PerspWindow(REAL(xlim), REAL(ylim), REAL(zlim), dd);
 
    /* Compute facet order. */
 
    PROTECT(depth = allocVector(REALSXP, (nrows(z) - 1)*(ncols(z) - 1)));
    PROTECT(indx = allocVector(INTSXP, (nrows(z) - 1)*(ncols(z) - 1)));
    DepthOrder(REAL(z), REAL(x), REAL(y), nrows(z), ncols(z),
	       REAL(depth), INTEGER(indx));
 
    /* Now we order the facets by depth and then draw them back to front.
     * This is the "painters" algorithm. */
 
    GMode(1, dd);
 
    if (dobox) {
        PerspBox(0, REAL(xlim), REAL(ylim), REAL(zlim), dd);
	if (doaxes) {
	    SEXP xl = STRING_ELT(xlab, 0), yl = STRING_ELT(ylab, 0),
		zl = STRING_ELT(zlab, 0);
	    PerspAxes(REAL(xlim), REAL(ylim), REAL(zlim),
		      (xl == NA_STRING)? "" : CHAR(xl), 
		      (yl == NA_STRING)? "" : CHAR(yl),
		      (zl == NA_STRING)? "" : CHAR(zl),
		      nTicks, tickType, dd);
	}
    }
 
    DrawFacets(REAL(z), REAL(x), REAL(y), nrows(z), ncols(z), INTEGER(indx),
	       1/xs, 1/ys, expand/zs,
	       INTEGER(col), ncol, INTEGER(border)[0]);
 
    if (dobox)
        PerspBox(1, REAL(xlim), REAL(ylim), REAL(zlim), dd);
    GMode(0, dd);
 
    GRestorePars(dd);
    UNPROTECT(10);
    if (GRecording(call))
        recordGraphicOperation(op, originalArgs, dd);
 
    PROTECT(x = allocVector(REALSXP, 16));
    PROTECT(y = allocVector(INTSXP, 2));
    for (i = 0; i < 4; i++)
	for (j = 0; j < 4; j++)
	    REAL(x)[i + j * 4] = VT[i][j];
    INTEGER(y)[0] = 4;
    INTEGER(y)[1] = 4;
    setAttrib(x, R_DimSymbol, y);
    UNPROTECT(2);
    return x;
}
 

Rev 29864	Rev 30691
1	`/*`	1	`/*`
2	`* R : A Computer Language for Statistical Data Analysis`	2	`* R : A Computer Language for Statistical Data Analysis`
3	`* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka`	3	`* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka`
4	`* Copyright (C) 1997-2002 Robert Gentleman, Ross Ihaka and the R core team.`	4	`* Copyright (C) 1997-2002 Robert Gentleman, Ross Ihaka and the R core team.`
5	`*`	5	`*`
6	`* This program is free software; you can redistribute it and/or modify`	6	`* This program is free software; you can redistribute it and/or modify`
7	`* it under the terms of the GNU General Public License as published by`	7	`* it under the terms of the GNU General Public License as published by`
8	`* the Free Software Foundation; either version 2 of the License, or`	8	`* the Free Software Foundation; either version 2 of the License, or`
9	`* (at your option) any later version.`	9	`* (at your option) any later version.`
10	`*`	10	`*`
11	`* This program is distributed in the hope that it will be useful,`	11	`* This program is distributed in the hope that it will be useful,`
12	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`	12	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
13	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`	13	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
14	`* GNU General Public License for more details.`	14	`* GNU General Public License for more details.`
15	`*`	15	`*`
16	`* You should have received a copy of the GNU General Public License`	16	`* You should have received a copy of the GNU General Public License`
17	`* along with this program; if not, write to the Free Software`	17	`* along with this program; if not, write to the Free Software`
18	`* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA`	18	`* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA`
19	`*/`	19	`*/`
20		20
21	`#ifdef HAVE_CONFIG_H`	21	`#ifdef HAVE_CONFIG_H`
22	`#include <config.h>`	22	`#include <config.h>`
23	`#endif`	23	`#endif`
24		24
25	`#include <Defn.h>`	25	`#include <Defn.h>`
26	`#include <Rmath.h>`	26	`#include <Rmath.h>`
27	`#include <Graphics.h>`	27	`#include <Graphics.h>`
28	`#include <Rdevices.h>`	28	`#include <Rdevices.h>`
29	`#include <Print.h>`	29	`#include <Print.h>`
30	`#include <R_ext/Boolean.h>`	30	`#include <R_ext/Boolean.h>`
31		31
32	`#ifndef HAVE_HYPOT`	32	`#ifndef HAVE_HYPOT`
33	`# define hypot pythag`	33	`# define hypot pythag`
34	`#endif`	34	`#endif`
35		35
36	`/* Conversion of degrees to radians */`	36	`/* Conversion of degrees to radians */`
37		37
38	`#define DegToRad(x) (0.01745329251994329576 * x)`	38	`#define DegToRad(x) (0.01745329251994329576 * x)`
39		39
40	`/* Definitions of data structures for vectors and */`	40	`/* Definitions of data structures for vectors and */`
41	`/* transformations in homogeneous 3d coordinates */`	41	`/* transformations in homogeneous 3d coordinates */`
42		42
43	`typedef double Vector3d[4];`	43	`typedef double Vector3d[4];`
44	`typedef double Trans3d[4][4];`	44	`typedef double Trans3d[4][4];`
45		45
46	`/* The viewing transformation matrix. */`	46	`/* The viewing transformation matrix. */`
47		47
48	`static SEXP gcall;`	48	`static SEXP gcall;`
49	`static Trans3d VT;`	49	`static Trans3d VT;`
50		50
51	`static void TransVector (Vector3d u, Trans3d T, Vector3d v)`	51	`static void TransVector (Vector3d u, Trans3d T, Vector3d v)`
52	`{`	52	`{`
53	`double sum;`	53	`double sum;`
54	`int i, j;`	54	`int i, j;`
55		55
56	`for (i = 0; i < 4; i++) {`	56	`for (i = 0; i < 4; i++) {`
57	`sum = 0;`	57	`sum = 0;`
58	`for (j = 0; j < 4; j++)`	58	`for (j = 0; j < 4; j++)`
59	`sum = sum + u[j] * T[j][i];`	59	`sum = sum + u[j] * T[j][i];`
60	`v[i] = sum;`	60	`v[i] = sum;`
61	`}`	61	`}`
62	`}`	62	`}`
63		63
64	`static void Accumulate (Trans3d T)`	64	`static void Accumulate (Trans3d T)`
65	`{`	65	`{`
66	`Trans3d U;`	66	`Trans3d U;`
67	`double sum;`	67	`double sum;`
68	`int i, j, k;`	68	`int i, j, k;`
69		69
70	`for (i = 0; i < 4; i++) {`	70	`for (i = 0; i < 4; i++) {`
71	`for (j = 0; j < 4; j++) {`	71	`for (j = 0; j < 4; j++) {`
72	`sum = 0;`	72	`sum = 0;`
73	`for (k = 0; k < 4; k++)`	73	`for (k = 0; k < 4; k++)`
74	`sum = sum + VT[i][k] * T[k][j];`	74	`sum = sum + VT[i][k] * T[k][j];`
75	`U[i][j] = sum;`	75	`U[i][j] = sum;`
76	`}`	76	`}`
77	`}`	77	`}`
78	`for (i = 0; i < 4; i++)`	78	`for (i = 0; i < 4; i++)`
79	`for (j = 0; j < 4; j++)`	79	`for (j = 0; j < 4; j++)`
80	`VT[i][j] = U[i][j];`	80	`VT[i][j] = U[i][j];`
81	`}`	81	`}`
82		82
83	`static void SetToIdentity (Trans3d T)`	83	`static void SetToIdentity (Trans3d T)`
84	`{`	84	`{`
85	`int i, j;`	85	`int i, j;`
86	`for (i = 0; i < 4; i++) {`	86	`for (i = 0; i < 4; i++) {`
87	`for (j = 0; j < 4; j++)`	87	`for (j = 0; j < 4; j++)`
88	`T[i][j] = 0;`	88	`T[i][j] = 0;`
89	`T[i][i] = 1;`	89	`T[i][i] = 1;`
90	`}`	90	`}`
91	`}`	91	`}`
92		92
93	`static void Translate (double x, double y, double z)`	93	`static void Translate (double x, double y, double z)`
94	`{`	94	`{`
95	`Trans3d T;`	95	`Trans3d T;`
96	`SetToIdentity(T);`	96	`SetToIdentity(T);`
97	`T[3][0] = x;`	97	`T[3][0] = x;`
98	`T[3][1] = y;`	98	`T[3][1] = y;`
99	`T[3][2] = z;`	99	`T[3][2] = z;`
100	`Accumulate(T);`	100	`Accumulate(T);`
101	`}`	101	`}`
102		102
103	`static void Scale (double x, double y, double z)`	103	`static void Scale (double x, double y, double z)`
104	`{`	104	`{`
105	`Trans3d T;`	105	`Trans3d T;`
106	`SetToIdentity(T);`	106	`SetToIdentity(T);`
107	`T[0][0] = x;`	107	`T[0][0] = x;`
108	`T[1][1] = y;`	108	`T[1][1] = y;`
109	`T[2][2] = z;`	109	`T[2][2] = z;`
110	`Accumulate(T);`	110	`Accumulate(T);`
111	`}`	111	`}`
112		112
113	`static void XRotate (double angle)`	113	`static void XRotate (double angle)`
114	`{`	114	`{`
115	`double c, s;`	115	`double c, s;`
116	`Trans3d T;`	116	`Trans3d T;`
117	`SetToIdentity(T);`	117	`SetToIdentity(T);`
118	`c = cos(DegToRad(angle));`	118	`c = cos(DegToRad(angle));`
119	`s = sin(DegToRad(angle));`	119	`s = sin(DegToRad(angle));`
120	`T[1][1] = c;`	120	`T[1][1] = c;`
121	`T[2][1] = -s;`	121	`T[2][1] = -s;`
122	`T[2][2] = c;`	122	`T[2][2] = c;`
123	`T[1][2] = s;`	123	`T[1][2] = s;`
124	`Accumulate(T);`	124	`Accumulate(T);`
125	`}`	125	`}`
126		126
127	`static void YRotate (double angle)`	127	`static void YRotate (double angle)`
128	`{`	128	`{`
129	`double c, s;`	129	`double c, s;`
130	`Trans3d T;`	130	`Trans3d T;`
131	`SetToIdentity(T);`	131	`SetToIdentity(T);`
132	`c = cos(DegToRad(angle));`	132	`c = cos(DegToRad(angle));`
133	`s = sin(DegToRad(angle));`	133	`s = sin(DegToRad(angle));`
134	`T[0][0] = c;`	134	`T[0][0] = c;`
135	`T[2][0] = s;`	135	`T[2][0] = s;`
136	`T[2][2] = c;`	136	`T[2][2] = c;`
137	`T[0][2] = -s;`	137	`T[0][2] = -s;`
138	`Accumulate(T);`	138	`Accumulate(T);`
139	`}`	139	`}`
140		140
141	`static void ZRotate (double angle)`	141	`static void ZRotate (double angle)`
142	`{`	142	`{`
143	`double c, s;`	143	`double c, s;`
144	`Trans3d T;`	144	`Trans3d T;`
145	`SetToIdentity(T);`	145	`SetToIdentity(T);`
146	`c = cos(DegToRad(angle));`	146	`c = cos(DegToRad(angle));`
147	`s = sin(DegToRad(angle));`	147	`s = sin(DegToRad(angle));`
148	`T[0][0] = c;`	148	`T[0][0] = c;`
149	`T[1][0] = -s;`	149	`T[1][0] = -s;`
150	`T[1][1] = c;`	150	`T[1][1] = c;`
151	`T[0][1] = s;`	151	`T[0][1] = s;`
152	`Accumulate(T);`	152	`Accumulate(T);`
153	`}`	153	`}`
154		154
155	`static void Perspective (double d)`	155	`static void Perspective (double d)`
156	`{`	156	`{`
157	`Trans3d T;`	157	`Trans3d T;`
158		158
159	`SetToIdentity(T);`	159	`SetToIdentity(T);`
160	`T[2][3] = -1 / d;`	160	`T[2][3] = -1 / d;`
161	`Accumulate(T);`	161	`Accumulate(T);`
162	`}`	162	`}`
163		163
164		164
165	`/* Stuff for labels on contour plots`	165	`/* Stuff for labels on contour plots`
166	`Originally written by Nicholas Hildreth`	166	`Originally written by Nicholas Hildreth`
167	`Adapted by Paul Murrell`	167	`Adapted by Paul Murrell`
168	`*/`	168	`*/`
169		169
170	`static`	170	`static`
171	`void FindCorners(double width, double height, SEXP label,`	171	`void FindCorners(double width, double height, SEXP label,`
172	`double x0, double y0, double x1, double y1,`	172	`double x0, double y0, double x1, double y1,`
173	`DevDesc *dd) {`	173	`DevDesc *dd) {`
174	`double delta = height / width;`	174	`double delta = height / width;`
175	`double dx = GConvertXUnits(x1 - x0, USER, INCHES, dd) * delta;`	175	`double dx = GConvertXUnits(x1 - x0, USER, INCHES, dd) * delta;`
176	`double dy = GConvertYUnits(y1 - y0, USER, INCHES, dd) * delta;`	176	`double dy = GConvertYUnits(y1 - y0, USER, INCHES, dd) * delta;`
177	`dx = GConvertYUnits(dx, INCHES, USER, dd);`	177	`dx = GConvertYUnits(dx, INCHES, USER, dd);`
178	`dy = GConvertXUnits(dy, INCHES, USER, dd);`	178	`dy = GConvertXUnits(dy, INCHES, USER, dd);`
179		179
180	`REAL(label)[0] = x0 + dy;`	180	`REAL(label)[0] = x0 + dy;`
181	`REAL(label)[4] = y0 - dx;`	181	`REAL(label)[4] = y0 - dx;`
182	`REAL(label)[1] = x0 - dy;`	182	`REAL(label)[1] = x0 - dy;`
183	`REAL(label)[5] = y0 + dx;`	183	`REAL(label)[5] = y0 + dx;`
184	`REAL(label)[3] = x1 + dy;`	184	`REAL(label)[3] = x1 + dy;`
185	`REAL(label)[7] = y1 - dx;`	185	`REAL(label)[7] = y1 - dx;`
186	`REAL(label)[2] = x1 - dy;`	186	`REAL(label)[2] = x1 - dy;`
187	`REAL(label)[6] = y1 + dx;`	187	`REAL(label)[6] = y1 + dx;`
188	`}`	188	`}`
189		189
190	`static`	190	`static`
191	`int TestLabelIntersection(SEXP label1, SEXP label2) {`	191	`int TestLabelIntersection(SEXP label1, SEXP label2) {`
192		192
193	`int i, j, l1, l2;`	193	`int i, j, l1, l2;`
194	`double Ax, Bx, Ay, By, ax, ay, bx, by, q1, q2;`	194	`double Ax, Bx, Ay, By, ax, ay, bx, by, q1, q2;`
195	`double dom;`	195	`double dom;`
196	`double result1, result2;`	196	`double result1, result2;`
197		197
198	`for (i = 0; i < 4; i++) {`	198	`for (i = 0; i < 4; i++) {`
199	`Ax = REAL(label1)[i];`	199	`Ax = REAL(label1)[i];`
200	`Ay = REAL(label1)[i+4];`	200	`Ay = REAL(label1)[i+4];`
201	`Bx = REAL(label1)[(i+1)%4];`	201	`Bx = REAL(label1)[(i+1)%4];`
202	`By = REAL(label1)[(i+1)%4+4];`	202	`By = REAL(label1)[(i+1)%4+4];`
203	`for (j = 0; j < 4; j++) {`	203	`for (j = 0; j < 4; j++) {`
204	`ax = REAL(label2)[j];`	204	`ax = REAL(label2)[j];`
205	`ay = REAL(label2)[j+4];`	205	`ay = REAL(label2)[j+4];`
206	`bx = REAL(label2)[(j+1)%4];`	206	`bx = REAL(label2)[(j+1)%4];`
207	`by = REAL(label2)[(j+1)%4+4];`	207	`by = REAL(label2)[(j+1)%4+4];`
208		208
209	`q1 = Ax*(ay-by);`	209	`q1 = Ax*(ay-by);`
210	`q2 = Ay*(bx-ax);`	210	`q2 = Ay*(bx-ax);`
211		211
212	`dom = Bxby - Bxay - Axby + Axay - bxBy + bxAy + axBy - axAy;`	212	`dom = Bxby - Bxay - Axby + Axay - bxBy + bxAy + axBy - axAy;`
213	`if (dom == 0.0) {`	213	`if (dom == 0.0) {`
214	`result1 = -1;`	214	`result1 = -1;`
215	`result2 = -1;`	215	`result2 = -1;`
216	`}`	216	`}`
217	`else {`	217	`else {`
218	`result1 = (bxAy - axAy - aybx - Axby + Axay + byax) / dom;`	218	`result1 = (bxAy - axAy - aybx - Axby + Axay + byax) / dom;`
219		219
220	`if (bx - ax == 0.0) {`	220	`if (bx - ax == 0.0) {`
221	`if (by - ay == 0.0)`	221	`if (by - ay == 0.0)`
222	`result2 = -1;`	222	`result2 = -1;`
223	`else`	223	`else`
224	`result2 = (Ay + (By - Ay) * result1 - ay) / (by - ay);`	224	`result2 = (Ay + (By - Ay) * result1 - ay) / (by - ay);`
225	`}`	225	`}`
226	`else`	226	`else`
227	`result2 = (Ax + (Bx - Ax) * result1 - ax) / (bx - ax);`	227	`result2 = (Ax + (Bx - Ax) * result1 - ax) / (bx - ax);`
228		228
229	`}`	229	`}`
230	`l1 = (result1 >= 0.0) && (result1 <= 1.0);`	230	`l1 = (result1 >= 0.0) && (result1 <= 1.0);`
231	`l2 = (result2 >= 0.0) && (result2 <= 1.0);`	231	`l2 = (result2 >= 0.0) && (result2 <= 1.0);`
232	`if (l1 && l2) return 1;`	232	`if (l1 && l2) return 1;`
233	`}`	233	`}`
234	`}`	234	`}`
235		235
236	`return 0;`	236	`return 0;`
237	`}`	237	`}`
238		238
239	`/* Checks whether a label window is inside view region */`	239	`/* Checks whether a label window is inside view region */`
240	`static int LabelInsideWindow(SEXP label, DevDesc *dd) {`	240	`static int LabelInsideWindow(SEXP label, DevDesc *dd) {`
241	`int i = 0;`	241	`int i = 0;`
242	`double x, y;`	242	`double x, y;`
243		243
244	`while (i < 4) {`	244	`while (i < 4) {`
245	`x = REAL(label)[i];`	245	`x = REAL(label)[i];`
246	`y = REAL(label)[i+4];`	246	`y = REAL(label)[i+4];`
247	`GConvert(&x, &y, USER, NDC, dd);`	247	`GConvert(&x, &y, USER, NDC, dd);`
248	`/* x = GConvertXUnits(REAL(label)[i], USER, NDC, dd);`	248	`/* x = GConvertXUnits(REAL(label)[i], USER, NDC, dd);`
249	`y = GConvertYUnits(REAL(label)[i+4], USER, NDC, dd); */`	249	`y = GConvertYUnits(REAL(label)[i+4], USER, NDC, dd); */`
250		250
251	`if ((x < 0) \|\| (x > 1) \|\|`	251	`if ((x < 0) \|\| (x > 1) \|\|`
252	`(y < 0) \|\| (y > 1))`	252	`(y < 0) \|\| (y > 1))`
253	`return 1;`	253	`return 1;`
254	`i += 1;`	254	`i += 1;`
255	`}`	255	`}`
256	`return 0;`	256	`return 0;`
257	`}`	257	`}`
258		258
259		259
260	`/* C o n t o u r P l o t t i n g */`	260	`/* C o n t o u r P l o t t i n g */`
261		261
262	`typedef struct SEG {`	262	`typedef struct SEG {`
263	`struct SEG *next;`	263	`struct SEG *next;`
264	`double x0;`	264	`double x0;`
265	`double y0;`	265	`double y0;`
266	`double x1;`	266	`double x1;`
267	`double y1;`	267	`double y1;`
268	`} SEG, *SEGP;`	268	`} SEG, *SEGP;`
269		269
270	`static SEGP *ctr_SegDB;`	270	`static SEGP *ctr_SegDB;`
271		271
272	`static int ctr_intersect(double z0, double z1, double zc, double *f)`	272	`static int ctr_intersect(double z0, double z1, double zc, double *f)`
273	`{`	273	`{`
274	`if ((z0 - zc) * (z1 - zc) < 0.0) {`	274	`if ((z0 - zc) * (z1 - zc) < 0.0) {`
275	`*f = (zc - z0) / (z1 - z0);`	275	`*f = (zc - z0) / (z1 - z0);`
276	`return 1;`	276	`return 1;`
277	`}`	277	`}`
278	`return 0;`	278	`return 0;`
279	`}`	279	`}`
280		280
281	`static SEGP ctr_newseg(double x0, double y0, double x1, double y1, SEGP prev)`	281	`static SEGP ctr_newseg(double x0, double y0, double x1, double y1, SEGP prev)`
282	`{`	282	`{`
283	`SEGP seg = (SEGP)R_alloc(1, sizeof(SEG));`	283	`SEGP seg = (SEGP)R_alloc(1, sizeof(SEG));`
284	`seg->x0 = x0;`	284	`seg->x0 = x0;`
285	`seg->y0 = y0;`	285	`seg->y0 = y0;`
286	`seg->x1 = x1;`	286	`seg->x1 = x1;`
287	`seg->y1 = y1;`	287	`seg->y1 = y1;`
288	`seg->next = prev;`	288	`seg->next = prev;`
289	`return seg;`	289	`return seg;`
290	`}`	290	`}`
291		291
292	`static void ctr_swapseg(SEGP seg)`	292	`static void ctr_swapseg(SEGP seg)`
293	`{`	293	`{`
294	`double x, y;`	294	`double x, y;`
295	`x = seg->x0;`	295	`x = seg->x0;`
296	`y = seg->y0;`	296	`y = seg->y0;`
297	`seg->x0 = seg->x1;`	297	`seg->x0 = seg->x1;`
298	`seg->y0 = seg->y1;`	298	`seg->y0 = seg->y1;`
299	`seg->x1 = x;`	299	`seg->x1 = x;`
300	`seg->y1 = y;`	300	`seg->y1 = y;`
301	`}`	301	`}`
302		302
303	`/* ctr_segdir(): Determine the entry direction to the next cell */`	303	`/* ctr_segdir(): Determine the entry direction to the next cell */`
304	`/* and update the cell indices */`	304	`/* and update the cell indices */`
305		305
306	`#define XMATCH(x0,x1) (fabs(x0-x1) == 0)`	306	`#define XMATCH(x0,x1) (fabs(x0-x1) == 0)`
307	`#define YMATCH(y0,y1) (fabs(y0-y1) == 0)`	307	`#define YMATCH(y0,y1) (fabs(y0-y1) == 0)`
308		308
309	`static int ctr_segdir(double xend, double yend, double x, double y,`	309	`static int ctr_segdir(double xend, double yend, double x, double y,`
310	`int i, int j, int nx, int ny)`	310	`int i, int j, int nx, int ny)`
311	`{`	311	`{`
312	`if (YMATCH(yend, y[*j])) {`	312	`if (YMATCH(yend, y[*j])) {`
313	`if (*j == 0)`	313	`if (*j == 0)`
314	`return 0;`	314	`return 0;`
315	`j = j - 1;`	315	`j = j - 1;`
316	`return 3;`	316	`return 3;`
317	`}`	317	`}`
318	`if (XMATCH(xend, x[*i])) {`	318	`if (XMATCH(xend, x[*i])) {`
319	`if (*i == 0)`	319	`if (*i == 0)`
320	`return 0;`	320	`return 0;`
321	`i = i - 1;`	321	`i = i - 1;`
322	`return 4;`	322	`return 4;`
323	`}`	323	`}`
324	`if (YMATCH(yend, y[*j + 1])) {`	324	`if (YMATCH(yend, y[*j + 1])) {`
325	`if (*j >= ny - 1)`	325	`if (*j >= ny - 1)`
326	`return 0;`	326	`return 0;`
327	`j = j + 1;`	327	`j = j + 1;`
328	`return 1;`	328	`return 1;`
329	`}`	329	`}`
330	`if (XMATCH(xend, x[*i + 1])) {`	330	`if (XMATCH(xend, x[*i + 1])) {`
331	`if (*i >= nx - 1)`	331	`if (*i >= nx - 1)`
332	`return 0;`	332	`return 0;`
333	`i = i + 1;`	333	`i = i + 1;`
334	`return 2;`	334	`return 2;`
335	`}`	335	`}`
336	`return 0;`	336	`return 0;`
337	`}`	337	`}`
338		338
339	`/* Search seglist for a segment with endpoint (xend, yend). */`	339	`/* Search seglist for a segment with endpoint (xend, yend). */`
340	`/* The cell entry direction is dir, and if tail=1/0 we are */`	340	`/* The cell entry direction is dir, and if tail=1/0 we are */`
341	`/* building the tail/head of a contour. The matching segment */`	341	`/* building the tail/head of a contour. The matching segment */`
342	`/* is pointed to by seg and the updated segment list (with */`	342	`/* is pointed to by seg and the updated segment list (with */`
343	`/* the matched segment stripped) is returned by the funtion. */`	343	`/* the matched segment stripped) is returned by the funtion. */`
344		344
345	`static SEGP ctr_segupdate(double xend, double yend, int dir, Rboolean tail,`	345	`static SEGP ctr_segupdate(double xend, double yend, int dir, Rboolean tail,`
346	`SEGP seglist, SEGP* seg)`	346	`SEGP seglist, SEGP* seg)`
347	`{`	347	`{`
348	`if (seglist == NULL) {`	348	`if (seglist == NULL) {`
349	`*seg = NULL;`	349	`*seg = NULL;`
350	`return NULL;`	350	`return NULL;`
351	`}`	351	`}`
352	`switch (dir) {`	352	`switch (dir) {`
353	`case 1:`	353	`case 1:`
354	`case 3:`	354	`case 3:`
355	`if (YMATCH(yend,seglist->y0)) {`	355	`if (YMATCH(yend,seglist->y0)) {`
356	`if (!tail)`	356	`if (!tail)`
357	`ctr_swapseg(seglist);`	357	`ctr_swapseg(seglist);`
358	`*seg = seglist;`	358	`*seg = seglist;`
359	`return seglist->next;`	359	`return seglist->next;`
360	`}`	360	`}`
361	`if (YMATCH(yend,seglist->y1)) {`	361	`if (YMATCH(yend,seglist->y1)) {`
362	`if (tail)`	362	`if (tail)`
363	`ctr_swapseg(seglist);`	363	`ctr_swapseg(seglist);`
364	`*seg = seglist;`	364	`*seg = seglist;`
365	`return seglist->next;`	365	`return seglist->next;`
366	`}`	366	`}`
367	`break;`	367	`break;`
368	`case 2:`	368	`case 2:`
369	`case 4:`	369	`case 4:`
370	`if (XMATCH(xend,seglist->x0)) {`	370	`if (XMATCH(xend,seglist->x0)) {`
371	`if (!tail)`	371	`if (!tail)`
372	`ctr_swapseg(seglist);`	372	`ctr_swapseg(seglist);`
373	`*seg = seglist;`	373	`*seg = seglist;`
374	`return seglist->next;`	374	`return seglist->next;`
375	`}`	375	`}`
376	`if (XMATCH(xend,seglist->x1)) {`	376	`if (XMATCH(xend,seglist->x1)) {`
377	`if (tail)`	377	`if (tail)`
378	`ctr_swapseg(seglist);`	378	`ctr_swapseg(seglist);`
379	`*seg = seglist;`	379	`*seg = seglist;`
380	`return seglist->next;`	380	`return seglist->next;`
381	`}`	381	`}`
382	`break;`	382	`break;`
383	`}`	383	`}`
384	`seglist->next = ctr_segupdate(xend, yend, dir, tail, seglist->next, seg);`	384	`seglist->next = ctr_segupdate(xend, yend, dir, tail, seglist->next, seg);`
385	`return seglist;`	385	`return seglist;`
386	`}`	386	`}`
387		387
388	`/* labelList, label1, and label2 are all SEXPs rather than being allocated`	388	`/* labelList, label1, and label2 are all SEXPs rather than being allocated`
389	`using R_alloc because they need to persist across calls to contour().`	389	`using R_alloc because they need to persist across calls to contour().`
390	`In do_contour() there is a vmaxget() ... vmaxset() around each call to`	390	`In do_contour() there is a vmaxget() ... vmaxset() around each call to`
391	`contour() to release all of the memory used in the drawing of the`	391	`contour() to release all of the memory used in the drawing of the`
392	`contour _lines_ at each contour level. We need to keep track of the`	392	`contour _lines_ at each contour level. We need to keep track of the`
393	`contour _labels_ for _all_ contour levels, hence we have to use a`	393	`contour _labels_ for _all_ contour levels, hence we have to use a`
394	`different memory allocation mechanism.`	394	`different memory allocation mechanism.`
395	`*/`	395	`*/`
396	`static SEXP labelList;`	396	`static SEXP labelList;`
397		397
398	`static`	398	`static`
399	`double distFromEdge(double xxx, double yyy, int iii, DevDesc *dd) {`	399	`double distFromEdge(double xxx, double yyy, int iii, DevDesc *dd) {`
400	`return fmin2(fmin2(xxx[iii]-Rf_gpptr(dd)->usr[0], Rf_gpptr(dd)->usr[1]-xxx[iii]),`	400	`return fmin2(fmin2(xxx[iii]-Rf_gpptr(dd)->usr[0], Rf_gpptr(dd)->usr[1]-xxx[iii]),`
401	`fmin2(yyy[iii]-Rf_gpptr(dd)->usr[2], Rf_gpptr(dd)->usr[3]-yyy[iii]));`	401	`fmin2(yyy[iii]-Rf_gpptr(dd)->usr[2], Rf_gpptr(dd)->usr[3]-yyy[iii]));`
402	`}`	402	`}`
403		403
404	`static`	404	`static`
405	`Rboolean useStart(double xxx, double yyy, int ns, DevDesc *dd) {`	405	`Rboolean useStart(double xxx, double yyy, int ns, DevDesc *dd) {`
406	`if (distFromEdge(xxx, yyy, 0, dd) < distFromEdge(xxx, yyy, ns-1, dd))`	406	`if (distFromEdge(xxx, yyy, 0, dd) < distFromEdge(xxx, yyy, ns-1, dd))`
407	`return TRUE;`	407	`return TRUE;`
408	`else`	408	`else`
409	`return FALSE;`	409	`return FALSE;`
410	`}`	410	`}`
411		411
412	`static`	412	`static`
413	`int findGapUp(double xxx, double yyy, int ns, double labelDistance,`	413	`int findGapUp(double xxx, double yyy, int ns, double labelDistance,`
414	`DevDesc *dd) {`	414	`DevDesc *dd) {`
415	`double dX, dY;`	415	`double dX, dY;`
416	`double dXC, dYC;`	416	`double dXC, dYC;`
417	`double distanceSum = 0;`	417	`double distanceSum = 0;`
418	`int n = 0;`	418	`int n = 0;`
419	`int jjj = 1;`	419	`int jjj = 1;`
420	`while ((jjj < ns) && (distanceSum < labelDistance)) {`	420	`while ((jjj < ns) && (distanceSum < labelDistance)) {`
421	`/* Find a gap big enough for the label`	421	`/* Find a gap big enough for the label`
422	`use several segments if necessary`	422	`use several segments if necessary`
423	`*/`	423	`*/`
424	`dX = xxx[jjj] - xxx[jjj - n - 1];`	424	`dX = xxx[jjj] - xxx[jjj - n - 1]; /* jjj - n - 1 == 0 */`
425	`dY = yyy[jjj] - yyy[jjj - n - 1];`	425	`dY = yyy[jjj] - yyy[jjj - n - 1];`
426	`dXC = GConvertXUnits(dX, USER, INCHES, dd);`	426	`dXC = GConvertXUnits(dX, USER, INCHES, dd);`
427	`dYC = GConvertYUnits(dY, USER, INCHES, dd);`	427	`dYC = GConvertYUnits(dY, USER, INCHES, dd);`
428	`distanceSum = hypot(dXC, dYC);`	428	`distanceSum = hypot(dXC, dYC);`
429	`jjj = (jjj + 1);`	429	`jjj++;`
430	`n += 1;`	430	`n++;`
431	`}`	431	`}`
432	`if (distanceSum < labelDistance)`	432	`if (distanceSum < labelDistance)`
433	`return 0;`	433	`return 0;`
434	`else`	434	`else`
435	`return n;`	435	`return n;`
436	`}`	436	`}`
437		437
438	`static`	438	`static`
439	`int findGapDown(double xxx, double yyy, int ns, double labelDistance,`	439	`int findGapDown(double xxx, double yyy, int ns, double labelDistance,`
440	`DevDesc *dd) {`	440	`DevDesc *dd) {`
441	`double dX, dY;`	441	`double dX, dY;`
442	`double dXC, dYC;`	442	`double dXC, dYC;`
443	`double distanceSum = 0;`	443	`double distanceSum = 0;`
444	`int n = 0;`	444	`int n = 0;`
445	`int jjj = ns - 2;`	445	`int jjj = ns - 2;`
446	`while ((jjj > -1) && (distanceSum < labelDistance)) {`	446	`while ((jjj > -1) && (distanceSum < labelDistance)) {`
447	`/* Find a gap big enough for the label`	447	`/* Find a gap big enough for the label`
448	`use several segments if necessary`	448	`use several segments if necessary`
449	`*/`	449	`*/`
450	`dX = xxx[jjj] - xxx[jjj + n + 1];`	450	`dX = xxx[jjj] - xxx[jjj + n + 1]; /jjj + n + 1 == ns -1 /`
451	`dY = yyy[jjj] - yyy[jjj + n + 1];`	451	`dY = yyy[jjj] - yyy[jjj + n + 1];`
452	`dXC = GConvertXUnits(dX, USER, INCHES, dd);`	452	`dXC = GConvertXUnits(dX, USER, INCHES, dd);`
453	`dYC = GConvertYUnits(dY, USER, INCHES, dd);`	453	`dYC = GConvertYUnits(dY, USER, INCHES, dd);`
454	`distanceSum = hypot(dXC, dYC);`	454	`distanceSum = hypot(dXC, dYC);`
455	`jjj--;`	455	`jjj--;`
456	`n -= 1;`	456	`n++;`
457	`}`	457	`}`
458	`if (distanceSum < labelDistance)`	458	`if (distanceSum < labelDistance)`
459	`return 0;`	459	`return 0;`
460	`else`	460	`else`
461	`return -n;`	461	`return n;`
462	`}`	462	`}`
463		463
464	`/*`	464	`/*`
465	`* Generate a list of segments for a single level`	465	`* Generate a list of segments for a single level`
466	`*/`	466	`*/`
467	`static SEGP* contourLines(double x, int nx, double y, int ny,`	467	`static SEGP* contourLines(double x, int nx, double y, int ny,`
468	`double *z, double zc, double atom)`	468	`double *z, double zc, double atom)`
469	`{`	469	`{`
470	`double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];`	470	`double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];`
471	`int i, j, k, l, m, nacode;`	471	`int i, j, k, l, m, nacode;`
472	`SEGP seglist;`	472	`SEGP seglist;`
473	`SEGP *segmentDB;`	473	`SEGP *segmentDB;`
474	`/* Initialize the segment data base */`	474	`/* Initialize the segment data base */`
475	`/* Note we must be careful about resetting */`	475	`/* Note we must be careful about resetting */`
476	`/* the top of the stack, otherwise we run out of */`	476	`/* the top of the stack, otherwise we run out of */`
477	`/* memory after a sequence of displaylist replays */`	477	`/* memory after a sequence of displaylist replays */`
478	`/*`	478	`/*`
479	`* This reset is done out in GEcontourLines`	479	`* This reset is done out in GEcontourLines`
480	`*/`	480	`*/`
481	`segmentDB = (SEGP)R_alloc(nxny, sizeof(SEGP));`	481	`segmentDB = (SEGP)R_alloc(nxny, sizeof(SEGP));`
482	`for (i = 0; i < nx; i++)`	482	`for (i = 0; i < nx; i++)`
483	`for (j = 0; j < ny; j++)`	483	`for (j = 0; j < ny; j++)`
484	`segmentDB[i + j * nx] = NULL;`	484	`segmentDB[i + j * nx] = NULL;`
485	`for (i = 0; i < nx - 1; i++) {`	485	`for (i = 0; i < nx - 1; i++) {`
486	`xl = x[i];`	486	`xl = x[i];`
487	`xh = x[i + 1];`	487	`xh = x[i + 1];`
488	`for (j = 0; j < ny - 1; j++) {`	488	`for (j = 0; j < ny - 1; j++) {`
489	`yl = y[j];`	489	`yl = y[j];`
490	`yh = y[j + 1];`	490	`yh = y[j + 1];`
491	`k = i + j * nx;`	491	`k = i + j * nx;`
492	`zll = z[k];`	492	`zll = z[k];`
493	`zhl = z[k + 1];`	493	`zhl = z[k + 1];`
494	`zlh = z[k + nx];`	494	`zlh = z[k + nx];`
495	`zhh = z[k + nx + 1];`	495	`zhh = z[k + nx + 1];`
496		496
497	`/* If the value at a corner is exactly equal to a contour level,`	497	`/* If the value at a corner is exactly equal to a contour level,`
498	`* change that value by a tiny amount */`	498	`* change that value by a tiny amount */`
499		499
500	`if (zll == zc) zll += atom;`	500	`if (zll == zc) zll += atom;`

The R Project SVN R

(root)/branches/djm-source/src/main/plot3d.c @ 34412 – Rev 29864 → 30691