The R Project SVN
R
(root)/trunk/src/main/plot3d.c – Rev 57770
Rev 55883 | Rev 59039 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
/** R : A Computer Language for Statistical Data Analysis* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka* Copyright (C) 1997--2011 The R Development 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, a copy is available at* http://www.r-project.org/Licenses/*/#ifdef HAVE_CONFIG_H#include <config.h>#endif#include <Defn.h>#include <float.h> /* for DBL_MAX */#include <Rmath.h>#include <Graphics.h>#include <Colors.h> /* for isNAcol */#include <Print.h>#include <R_ext/Boolean.h>/* Conversion of degrees to radians */#define DegToRad(x) (DEG2RAD * 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 plotsOriginally written by Nicholas HildrethAdapted by Paul Murrell*/staticvoid FindCorners(double width, double height, SEXP label,double x0, double y0, double x1, double y1,pGEDevDesc 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;}staticint TestLabelIntersection(SEXP label1, SEXP label2) {int i, j, l1, l2;double Ax, Bx, Ay, By, ax, ay, bx, by;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];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;elseresult2 = (Ay + (By - Ay) * result1 - ay) / (by - ay);}elseresult2 = (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, pGEDevDesc 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 allocatedusing R_alloc because they need to persist across calls to contour().In do_contour() there is a vmaxget() ... vmaxset() around each call tocontour() to release all of the memory used in the drawing of thecontour _lines_ at each contour level. We need to keep track of thecontour _labels_ for _all_ contour levels, hence we have to use adifferent memory allocation mechanism.*/static SEXP labelList;staticdouble distFromEdge(double *xxx, double *yyy, int iii, pGEDevDesc dd) {return fmin2(fmin2(xxx[iii]-gpptr(dd)->usr[0], gpptr(dd)->usr[1]-xxx[iii]),fmin2(yyy[iii]-gpptr(dd)->usr[2], gpptr(dd)->usr[3]-yyy[iii]));}staticRboolean useStart(double *xxx, double *yyy, int ns, pGEDevDesc dd) {if (distFromEdge(xxx, yyy, 0, dd) < distFromEdge(xxx, yyy, ns-1, dd))return TRUE;elsereturn FALSE;}staticint findGapUp(double *xxx, double *yyy, int ns, double labelDistance,pGEDevDesc 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 labeluse 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;elsereturn n;}staticint findGapDown(double *xxx, double *yyy, int ns, double labelDistance,pGEDevDesc 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 labeluse 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;elsereturn 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;}#define CONTOUR_LIST_STEP 100#define CONTOUR_LIST_LEVEL 0#define CONTOUR_LIST_X 1#define CONTOUR_LIST_Y 2static 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*/staticint 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, 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_contour_segments: prevent inf.loop (shouldn't be needed) */while (s && ns < max_contour_segments) {ns++;s = s->next;}if(ns == max_contour_segments)warning(_("contour(): circular/long seglist -- set %s > %d?"),"options(\"max.contour.segments\")", max_contour_segments);/** "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_contour_segments) {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.* Apart from here, used in package clines.*/SEXP GEcontourLines(double *x, int nx, double *y, int ny,double *z, double *levels, int nl){const void *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"));elsewarning(_("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_contour_segments 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 */} elsemainlist = VECTOR_ELT(container, 0);UNPROTECT(1); /* UNPROTECT container */return mainlist;}SEXP attribute_hidden do_contourLines(SEXP call, SEXP op, SEXP args, SEXP env){SEXP c, x, y, z;int nx, ny, nc;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);}/** 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,double atom, pGEDevDesc dd){/* draw a contour for one given contour level 'zc' */const void *vmax;double xend, yend;int i, ii, j, jj, ns, dir;SEGP seglist, seg, s, start, end;double *xxx, *yyy;double variance, dX, dY, deltaX, deltaY;double dXC, dYC;int range=0, indx=0, n; /* -Wall */double lowestVariance;double squareSum;int iii, jjj;double distanceSum, labelDistance, avgGradient;char buffer[255];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_contourRprintf("contour(lev = %g):\n", zc);#endifvmax = 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_contour_segments: prevent inf.loop (shouldn't be needed) */while (s && ns < max_contour_segments) {ns++;s = s->next;}if(ns == max_contour_segments)warning(_("contour(): circular/long seglist -- set %s > %d?"),"options(\"max.contour.segments\")", max_contour_segments);/* contour midpoint : use for labelling sometime (not yet!)int ns2;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_contour_segments) {s = s->next;xxx[ns] = s->x0;yyy[ns++] = s->y0;}xxx[ns] = s->x1;yyy[ns++] = s->y1;#ifdef DEBUG_contourRprintf(" [%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]);elsefor(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 themOtherwise stringify the z-value of the contour */cetype_t enc = CE_NATIVE;buffer[0] = ' ';if (!isNull(labels)) {int numl = length(labels);strcpy(&buffer[1], CHAR(STRING_ELT(labels, cnum % numl)));enc = getCharCE(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))); /* ASCII */UNPROTECT(1);}buffer[strlen(buffer)+1] = '\0';buffer[strlen(buffer)] = ' ';labelDistance = GStrWidth(buffer, enc, INCHES, dd);labelHeight = GStrHeight(buffer, enc, INCHES, dd);if (labelDistance > 0) {/* Try to find somewhere to draw the label */switch (method) {case 0: /* draw label at one end of contouroverwriting contour line*/if (useStart(xxx, yyy, ns, dd) )indx = 0;elseindx = ns - 1;break;case 1: /* draw label at one end of contourembedded in contourno 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 contourembedded in contour lineno 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;n = 0;jjj = (iii + 1);while ((jjj < ns-1) &&(distanceSum < labelDistance)) {/* Find a gap big enough for the labeluse 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 gradientsof the segments that will make way for thelabel*/deltaX = xxx[jjj] - xxx[jjj - 1];deltaY = yyy[jjj] - yyy[jjj - 1];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;}}if (lowestVariance < 9999999)gotLabel = TRUE;}} /* switch (method) */if (method == 0) {GPolyline(ns, xxx, yyy, USER, dd);GText(xxx[indx], yyy[indx], USER, buffer,CE_NATIVE/*FIX*/,.5, .5, 0, dd);}else {if (indx > 0)GPolyline(indx+1, xxx, yyy, USER, dd);if (ns-1-indx-range > 0)GPolyline(ns-indx-range, xxx+indx+range, yyy+indx+range,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] - gpptr(dd)->usr[0]),(gpptr(dd)->usr[1] - xxx[indx]));dx2 = fmin2((gpptr(dd)->usr[1] - xxx[indx+range]),(xxx[indx+range] - gpptr(dd)->usr[0]));if (dx1 < dx2) {closest = 0;dmin = dx1;} else {closest = 1;dmin = dx2;}dy1 = fmin2((yyy[indx] - gpptr(dd)->usr[2]),(gpptr(dd)->usr[3] - yyy[indx]));if (closest && (dy1 < dmin)) {closest = 0;dmin = dy1;} else if (dy1 < dmin)dmin = dy1;dy2 = fmin2((gpptr(dd)->usr[3] - yyy[indx+range]),(yyy[indx+range] - 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 labeldoesn'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 ofangle to draw text*/GConvert(&ux, &uy, USER, INCHES, dd);GConvert(&vx, &vy, USER, INCHES, dd);/* 0, .5 => left, centre justified */GText (ux, uy, INCHES, buffer,CE_NATIVE/*FIX*/,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 thingprotected is likely labelList, and that needsto be preserved across calls */}SEXP attribute_hidden Rg_contourDef(void){return ScalarLogical(GEcurrentDevice()->dev->useRotatedTextInContour);}/* contour(x, y, z, levels, labels, labcex, drawlabels,* method, vfont, col = col, lty = lty, lwd = lwd)*/SEXP attribute_hidden 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, lwdsave, fontsave = 1 /* -Wall */;rcolor colsave;double cexsave;double atom, zmin, zmax;const void *vmax, *vmax0;char familysave[201];int method;Rboolean drawLabels;double labcex;pGEDevDesc dd = GEcurrentDevice();SEXP result = R_NilValue;GCheckState(dd);if (length(args) < 4)error(_("too few arguments"));PrintDefaults(); /* prepare for labelformat */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)error(_("invalid '%s' value"), "method");PROTECT(vfont = FixupVFont(CAR(args)));if (!isNull(vfont)) {strncpy(familysave, gpptr(dd)->family, 201);strncpy(gpptr(dd)->family, "Her ", 201);gpptr(dd)->family[3] = INTEGER(vfont)[0];fontsave = gpptr(dd)->font;gpptr(dd)->font = INTEGER(vfont)[1];}args = CDR(args);rawcol = CAR(args);PROTECT(col = FixupCol(rawcol, R_TRANWHITE));ncol = length(col);args = CDR(args);PROTECT(lty = FixupLty(CAR(args), gpptr(dd)->lty));nlty = length(lty);args = CDR(args);PROTECT(lwd = FixupLwd(CAR(args), gpptr(dd)->lwd));nlwd = length(lwd);args = CDR(args);if (nx < 2 || ny < 2)error(_("insufficient 'x' or 'y' values"));if (nrows(z) != nx || ncols(z) != ny)error(_("dimension mismatch"));if (nc < 1)error(_("no contour values"));for (i = 0; i < nx; i++) {if (!R_FINITE(REAL(x)[i]))error(_("missing 'x' values"));if (i > 0 && REAL(x)[i] < REAL(x)[i - 1])error(_("increasing 'x' values expected"));}for (i = 0; i < ny; i++) {if (!R_FINITE(REAL(y)[i]))error(_("missing 'y' values"));if (i > 0 && REAL(y)[i] < REAL(y)[i - 1])error(_("increasing 'y' values expected"));}for (i = 0; i < nc; i++)if (!R_FINITE(REAL(c)[i]))error(_("invalid 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"));elsewarning(_("all z values are NA"));UNPROTECT(4);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_contour_segments 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 = gpptr(dd)->lty;colsave = gpptr(dd)->col;lwdsave = gpptr(dd)->lwd;cexsave = gpptr(dd)->cex;labelList = PROTECT(R_NilValue);/* draw contour for levels[i] */GMode(1, dd);for (i = 0; i < nc; i++) {vmax = vmaxget();gpptr(dd)->lty = INTEGER(lty)[i % nlty];if (gpptr(dd)->lty == NA_INTEGER)gpptr(dd)->lty = ltysave;if (isNAcol(rawcol, i, ncol))gpptr(dd)->col = colsave;elsegpptr(dd)->col = INTEGER(col)[i % ncol];gpptr(dd)->lwd = REAL(lwd)[i % nlwd];if (!R_FINITE(gpptr(dd)->lwd))gpptr(dd)->lwd = lwdsave;gpptr(dd)->cex = labcex;contour(x, nx, y, ny, z, REAL(c)[i], labels, i,drawLabels, method - 1, atom, dd);vmaxset(vmax);}GMode(0, dd);vmaxset(vmax0);gpptr(dd)->lty = ltysave;gpptr(dd)->col = colsave;gpptr(dd)->lwd = lwdsave;gpptr(dd)->cex = cexsave;if(!isNull(vfont)) {strncpy(gpptr(dd)->family, familysave, 201);gpptr(dd)->font = fontsave;}UNPROTECT(5);/* NOTE: only record operation if no "error" *//* NOTE: on replay, call == R_NilValue */if (GRecording(call, dd))GErecordGraphicOperation(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 voidFindCutPoints(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 voidFindPolygonVertices(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 ofdo_contour, do_filledcontour & do_image ...*//* filledcontour(x, y, z, levels, col) */SEXP attribute_hidden do_filledcontour(SEXP call, SEXP op, SEXP args, SEXP env){SEXP oargs, sx, sy, sz, sc, scol;double *x, *y, *z, *c;rcolor *col;int i, j, k, npt, nx, ny, nc, ncol, colsave, xpdsave;double px[8], py[8], pz[8];pGEDevDesc dd = GEcurrentDevice();GCheckState(dd);checkArity(op,args);PrintDefaults(); /* prepare for labelformat */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);args = CDR(args);sc = CAR(args);/* levels */internalTypeCheck(call, sc, REALSXP);nc = length(sc);args = CDR(args);if (nx < 2 || ny < 2)error(_("insufficient 'x' or 'y' values"));if (nrows(sz) != nx || ncols(sz) != ny)error(_("dimension mismatch"));if (nc < 1)error(_("no contour values"));PROTECT(scol = FixupCol(CAR(args), R_TRANWHITE));ncol = length(scol);/* Shorthand Pointers */x = REAL(sx);y = REAL(sy);z = REAL(sz);c = REAL(sc);col = (rcolor *) 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 = gpptr(dd)->col;xpdsave = gpptr(dd)->xpd;/* override par("xpd") and force clipping to plot region */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],R_TRANWHITE, dd);}}}GMode(0, dd);gpptr(dd)->col = colsave;gpptr(dd)->xpd = xpdsave;UNPROTECT(1);if (GRecording(call, dd))GErecordGraphicOperation(op, oargs, dd);return R_NilValue;badxy:error(_("invalid x / y values or limits"));badlev:error(_("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 attribute_hidden 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, xpdsave;rcolor colsave;pGEDevDesc dd = GEcurrentDevice();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), R_TRANWHITE));nc = LENGTH(sc);/* Shorthand Pointers */x = REAL(sx);y = REAL(sy);z = INTEGER(sz);c = (unsigned*)INTEGER(sc);/* Check of grid coordinates now done in C code */colsave = gpptr(dd)->col;xpdsave = gpptr(dd)->xpd;/* override par("xpd") and force clipping to plot region */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],R_TRANWHITE, dd);}}GMode(0, dd);gpptr(dd)->col = colsave;gpptr(dd)->xpd = xpdsave;UNPROTECT(1);if (GRecording(call, dd))GErecordGraphicOperation(op, oargs, dd);return R_NilValue;}/* 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 ensurethat 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;pGEDevDesc dd;dd = GEcurrentDevice();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);}}}static void PerspWindow(double *xlim, double *ylim, double *zlim, pGEDevDesc 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 teston the surfaces of the xlim/ylim/zlim box around the plot.If front = 0, only the faces with their inside toward theeyepoint are drawn. If front = 1, only the faces withtheir outside toward the eye are drawn. This lets us carryout hidden line removal by drawing any faces which will beobscured before the surface, and those which will not beobscured after the surface.Unfortunately as PR#202 showed, this is simplistic as the surfacecan go outside the box.*//* 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},};static void PerspBox(int front, double *x, double *y, double *z,char *EdgeDone, pGEDevDesc 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 = gpptr(dd)->lty;gpptr(dd)->lty = front ? LTY_DOTTED : 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);}}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;elsedy = y1 - y2;if (dx == 0) {if (dy > 0)angle = 90;elseangle = 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,const char *label, cetype_t enc, pGEDevDesc dd){Vector3d u1={0.,0.,0.,0.}, u2={0.,0.,0.,0.}, u3={0.,0.,0.,0.}, 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;double cexsave = gpptr(dd)->cex;int fontsave = gpptr(dd)->font;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" ticksin order to leave room for the tick labels */switch (tickType) {case 1: /* "simple": just an arrow parallel to axis, indicating directionof 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 *//* change in 2.5.0 to use cex.lab and font.lab */gpptr(dd)->cex = gpptr(dd)->cexbase * gpptr(dd)->cexlab;gpptr(dd)->font = gpptr(dd)->fontlab;GText(v3[0]/v3[3], v3[1]/v3[3], USER, label, enc, .5, .5,labelAngle(v1[0]/v1[3], v1[1]/v1[3], v2[0]/v2[3], v2[1]/v2[3]),dd);/* Draw axis ticks *//* change in 2.5.0 to use cex.axis and font.axis */gpptr(dd)->cex = gpptr(dd)->cexbase * gpptr(dd)->cexaxis;gpptr(dd)->font = gpptr(dd)->fontaxis;switch (tickType) {case 1: /* "simple": just an arrow parallel to axis, indicating directionof 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)),getCharCE(STRING_ELT(lab, i)),.5, .5, 0, dd);}UNPROTECT(2);break;}gpptr(dd)->cex = cexsave;gpptr(dd)->font = fontsave;}/* 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,const char *xlab, cetype_t xenc,const char *ylab, cetype_t yenc,const char *zlab, cetype_t zenc,int nTicks, int tickType, pGEDevDesc 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 = gpptr(dd)->xpd;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;} elsewarning(_("Axis orientation not calculated"));PerspAxis(x, y, z, xAxis, 0, nTicks, tickType, xlab, xenc, dd);PerspAxis(x, y, z, yAxis, 1, nTicks, tickType, ylab, yenc, 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;} elsewarning(_("Axis orientation not calculated"));PerspAxis(x, y, z, zAxis, 2, nTicks, tickType, zlab, zenc, dd);gpptr(dd)->xpd = xpdsave;}SEXP attribute_hidden 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 = 0.0, yc = 0.0, zc = 0.0, xs = 0.0, ys = 0.0, zs = 0.0;int i, j, scale, ncol, dobox, doaxes, nTicks, tickType;char EdgeDone[12]; /* Which edges have been drawn previously */pGEDevDesc dd;if (length(args) < 24) /* 24 plus any inline par()s */error(_("too few parameters"));gcall = call;originalArgs = args;PROTECT(x = coerceVector(CAR(args), REALSXP));if (length(x) < 2) error(_("invalid '%s' argument"), "x");args = CDR(args);PROTECT(y = coerceVector(CAR(args), REALSXP));if (length(y) < 2) error(_("invalid '%s' argument"), "y");args = CDR(args);PROTECT(z = coerceVector(CAR(args), REALSXP));if (!isMatrix(z) || nrows(z) != length(x) || ncols(z) != length(y))error(_("invalid '%s' argument"), "z");args = CDR(args);PROTECT(xlim = coerceVector(CAR(args), REALSXP));if (length(xlim) != 2) error(_("invalid '%s' argument"), "xlim");args = CDR(args);PROTECT(ylim = coerceVector(CAR(args), REALSXP));if (length(ylim) != 2) error(_("invalid '%s' argument"), "ylim");args = CDR(args);PROTECT(zlim = coerceVector(CAR(args), REALSXP));if (length(zlim) != 2) error(_("invalid '%s' argument"), "zlim");args = CDR(args);/* Checks on x/y/z Limits */if (!LimitCheck(REAL(xlim), &xc, &xs))error(_("invalid 'x' limits"));if (!LimitCheck(REAL(ylim), &yc, &ys))error(_("invalid 'y' limits"));if (!LimitCheck(REAL(zlim), &zc, &zs))error(_("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;elseDoLighting = 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)error(_("invalid viewing parameters"));if (!R_FINITE(expand) || expand < 0)error(_("invalid '%s' value"), "expand");if (scale == NA_LOGICAL)scale = 0;if ((nTicks == NA_INTEGER) || (nTicks < 0))error(_("invalid '%s' value"), "nticks");if ((tickType == NA_INTEGER) || (tickType < 1) || (tickType > 2))error(_("invalid '%s' value"), "ticktype");dd = GEcurrentDevice();GNewPlot(GRecording(call, dd));PROTECT(col = FixupCol(col, gpptr(dd)->bg));ncol = LENGTH(col);if (ncol < 1) error(_("invalid '%s' specification"), "col");if(!R_OPAQUE(INTEGER(col)[0])) DoLighting = FALSE;PROTECT(border = FixupCol(border, gpptr(dd)->fg));if (length(border) < 1)error(_("invalid '%s' specification"), "border");GSetState(1, dd);GSavePars(dd);ProcessInlinePars(args, dd, call);if (length(border) > 1)gpptr(dd)->fg = INTEGER(border)[0];gpptr(dd)->xlog = 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:We order the facets by depth and then draw them back to front.This is the "painters" algorithm. */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));GMode(1, dd);if (dobox) {/* Draw (solid) faces which face away from the viewer */PerspBox(0, REAL(xlim), REAL(ylim), REAL(zlim), EdgeDone, 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), getCharCE(xl),(yl == NA_STRING) ? "" : CHAR(yl), getCharCE(yl),(zl == NA_STRING) ? "" : CHAR(zl), getCharCE(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]);/* Draw (dotted) not-already-plotted edges of faces which facetowards from the viewer */if (dobox)PerspBox(1, REAL(xlim), REAL(ylim), REAL(zlim), EdgeDone, dd);GMode(0, dd);GRestorePars(dd);UNPROTECT(10);if (GRecording(call, dd))GErecordGraphicOperation(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;}