The R Project SVN R

Rev

Rev 27215 | Rev 29959 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1284 ihaka 1 
/*
2 
 *  R : A Computer Language for Statistical Data Analysis
3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
19134 ripley 4 
 *  Copyright (C) 1997-2002   Robert Gentleman, Ross Ihaka and the R core team.
1284 ihaka 5 
 *
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
8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
10 
 *
11 
 *  This program is distributed in the hope that it will be useful,
12 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 
 *  GNU General Public License for more details.
15 
 *
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
5458 ripley 18 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
1284 ihaka 19 
 */
20 
 
5187 hornik 21 
#ifdef HAVE_CONFIG_H
7701 hornik 22 
#include <config.h>
5187 hornik 23 
#endif
24 
 
11499 ripley 25 
#include <Defn.h>
26 
#include <Rmath.h>
27 
#include <Graphics.h>
12778 pd 28 
#include <Rdevices.h>
11499 ripley 29 
#include <Print.h>
13792 pd 30 
#include <R_ext/Boolean.h>
1284 ihaka 31 
 
9417 ripley 32 
#ifndef HAVE_HYPOT
33 
# define hypot pythag
34 
#endif
7813 murrell 35 
 
36 
/* Conversion of degrees to radians */
37 
 
38 
#define DegToRad(x) (0.01745329251994329576 * x)
39 
 
40 
/* Definitions of data structures for vectors and */
41 
/* transformations in homogeneous 3d coordinates */
42 
 
43 
typedef double Vector3d[4];
44 
typedef double Trans3d[4][4];
45 
 
46 
/* The viewing transformation matrix. */
47 
 
48 
static SEXP gcall;
49 
static Trans3d VT;
50 
 
51 
static void TransVector (Vector3d u, Trans3d T, Vector3d v)
52 
{
53 
    double sum;
54 
    int i, j;
55 
 
56 
    for (i = 0; i < 4; i++) {
57 
	sum = 0;
58 
	for (j = 0; j < 4; j++)
59 
	    sum = sum + u[j] * T[j][i];
60 
	v[i] = sum;
61 
    }
62 
}
63 
 
64 
static void Accumulate (Trans3d T)
65 
{
66 
    Trans3d U;
67 
    double sum;
68 
    int i, j, k;
69 
 
70 
    for (i = 0; i < 4; i++) {
71 
	for (j = 0; j < 4; j++) {
72 
	    sum = 0;
73 
	    for (k = 0; k < 4; k++)
74 
		sum = sum + VT[i][k] * T[k][j];
75 
	    U[i][j] = sum;
76 
	}
77 
    }
78 
    for (i = 0; i < 4; i++)
79 
	for (j = 0; j < 4; j++)
80 
	    VT[i][j] = U[i][j];
81 
}
82 
 
83 
static void SetToIdentity (Trans3d T)
84 
{
85 
    int i, j;
86 
    for (i = 0; i < 4; i++) {
87 
	for (j = 0; j < 4; j++)
88 
	    T[i][j] = 0;
89 
	T[i][i] = 1;
90 
    }
91 
}
92 
 
93 
static void Translate (double x, double y, double z)
94 
{
95 
    Trans3d T;
96 
    SetToIdentity(T);
97 
    T[3][0] = x;
98 
    T[3][1] = y;
99 
    T[3][2] = z;
100 
    Accumulate(T);
101 
}
102 
 
103 
static void Scale (double x, double y, double z)
104 
{
105 
    Trans3d T;
106 
    SetToIdentity(T);
107 
    T[0][0] = x;
108 
    T[1][1] = y;
109 
    T[2][2] = z;
110 
    Accumulate(T);
111 
}
112 
 
113 
static void XRotate (double angle)
114 
{
115 
    double c, s;
116 
    Trans3d T;
117 
    SetToIdentity(T);
118 
    c = cos(DegToRad(angle));
119 
    s = sin(DegToRad(angle));
120 
    T[1][1] = c;
121 
    T[2][1] = -s;
122 
    T[2][2] = c;
123 
    T[1][2] = s;
124 
    Accumulate(T);
125 
}
126 
 
127 
static void YRotate (double angle)
128 
{
129 
    double c, s;
130 
    Trans3d T;
131 
    SetToIdentity(T);
132 
    c = cos(DegToRad(angle));
133 
    s = sin(DegToRad(angle));
134 
    T[0][0] = c;
135 
    T[2][0] = s;
136 
    T[2][2] = c;
137 
    T[0][2] = -s;
138 
    Accumulate(T);
139 
}
140 
 
141 
static void ZRotate (double angle)
142 
{
143 
    double c, s;
144 
    Trans3d T;
145 
    SetToIdentity(T);
146 
    c = cos(DegToRad(angle));
147 
    s = sin(DegToRad(angle));
148 
    T[0][0] = c;
149 
    T[1][0] = -s;
150 
    T[1][1] = c;
151 
    T[0][1] = s;
152 
    Accumulate(T);
153 
}
154 
 
155 
static void Perspective (double d)
156 
{
157 
    Trans3d T;
158 
 
159 
    SetToIdentity(T);
160 
    T[2][3] = -1 / d;
161 
    Accumulate(T);
162 
}
163 
 
164 
 
9279 maechler 165 
/* Stuff for labels on contour plots
7813 murrell 166 
   Originally written by Nicholas Hildreth
167 
   Adapted by Paul Murrell
168 
*/
169 
 
7824 ripley 170 
static
7813 murrell 171 
void FindCorners(double width, double height, SEXP label,
172 
		 double x0, double y0, double x1, double y1,
173 
		 DevDesc *dd) {
174 
    double delta = height / width;
175 
    double dx = GConvertXUnits(x1 - x0, USER, INCHES, dd) * delta;
176 
    double dy = GConvertYUnits(y1 - y0, USER, INCHES, dd) * delta;
177 
    dx = GConvertYUnits(dx, INCHES, USER, dd);
178 
    dy = GConvertXUnits(dy, INCHES, USER, dd);
9279 maechler 179 
 
7813 murrell 180 
    REAL(label)[0] = x0 + dy;
181 
    REAL(label)[4] = y0 - dx;
182 
    REAL(label)[1] = x0 - dy;
183 
    REAL(label)[5] = y0 + dx;
184 
    REAL(label)[3] = x1 + dy;
185 
    REAL(label)[7] = y1 - dx;
186 
    REAL(label)[2] = x1 - dy;
9279 maechler 187 
    REAL(label)[6] = y1 + dx;
7813 murrell 188 
}
189 
 
7824 ripley 190 
static
7813 murrell 191 
int TestLabelIntersection(SEXP label1, SEXP label2) {
9279 maechler 192 
 
9487 ripley 193 
    int i, j, l1, l2;
7813 murrell 194 
    double Ax, Bx, Ay, By, ax, ay, bx, by, q1, q2;
7824 ripley 195 
    double dom;
7813 murrell 196 
    double result1, result2;
10886 maechler 197 
 
9487 ripley 198 
    for (i = 0; i < 4; i++) {
7813 murrell 199 
	Ax = REAL(label1)[i];
200 
	Ay = REAL(label1)[i+4];
201 
	Bx = REAL(label1)[(i+1)%4];
202 
	By = REAL(label1)[(i+1)%4+4];
9487 ripley 203 
	for (j = 0; j < 4; j++) {
7813 murrell 204 
	    ax = REAL(label2)[j];
205 
	    ay = REAL(label2)[j+4];
206 
	    bx = REAL(label2)[(j+1)%4];
207 
	    by = REAL(label2)[(j+1)%4+4];
9279 maechler 208 
 
7813 murrell 209 
	    q1 = Ax*(ay-by);
210 
	    q2 = Ay*(bx-ax);
9279 maechler 211 
 
7813 murrell 212 
	    dom = Bx*by - Bx*ay - Ax*by + Ax*ay - bx*By + bx*Ay + ax*By - ax*Ay;
213 
	    if (dom == 0.0) {
214 
		result1 = -1;
215 
		result2 = -1;
216 
	    }
217 
	    else {
218 
		result1 = (bx*Ay - ax*Ay - ay*bx - Ax*by + Ax*ay + by*ax) / dom;
9279 maechler 219 
 
7813 murrell 220 
		if (bx - ax == 0.0) {
221 
		    if (by - ay == 0.0)
222 
			result2 = -1;
223 
		    else
224 
			result2 = (Ay + (By - Ay) * result1 - ay) / (by - ay);
225 
		}
226 
		else
227 
		    result2 = (Ax + (Bx - Ax) * result1 - ax) / (bx - ax);
9279 maechler 228 
 
7813 murrell 229 
	    }
9487 ripley 230 
	    l1 = (result1 >= 0.0) && (result1 <= 1.0);
231 
	    l2 = (result2 >= 0.0) && (result2 <= 1.0);
232 
	    if (l1 && l2) return 1;
7813 murrell 233 
	}
234 
    }
9279 maechler 235 
 
7813 murrell 236 
    return 0;
237 
}
238 
 
239 
/*** Checks whether a label window is inside view region ***/
7824 ripley 240 
static int LabelInsideWindow(SEXP label, DevDesc *dd) {
7813 murrell 241 
    int i = 0;
242 
    double x, y;
9279 maechler 243 
 
7813 murrell 244 
    while (i < 4) {
245 
	x = REAL(label)[i];
246 
	y = REAL(label)[i+4];
247 
	GConvert(&x, &y, USER, NDC, dd);
8062 murrell 248 
	/*	x = GConvertXUnits(REAL(label)[i], USER, NDC, dd);
249 
		y = GConvertYUnits(REAL(label)[i+4], USER, NDC, dd); */
9279 maechler 250 
 
7813 murrell 251 
	if ((x < 0) || (x > 1) ||
252 
	    (y < 0) || (y > 1))
253 
	    return 1;
254 
	i += 1;
255 
    }
256 
    return 0;
257 
}
258 
 
259 
 
1284 ihaka 260 
	/*  C o n t o u r   P l o t t i n g  */
261 
 
262 
typedef struct SEG {
263 
    struct SEG *next;
264 
    double x0;
265 
    double y0;
266 
    double x1;
267 
    double y1;
268 
} SEG, *SEGP;
269 
 
270 
static SEGP *ctr_SegDB;
271 
 
272 
static int ctr_intersect(double z0, double z1, double zc, double *f)
273 
{
1839 ihaka 274 
    if ((z0 - zc) * (z1 - zc) < 0.0) {
1284 ihaka 275 
	*f = (zc - z0) / (z1 -	z0);
276 
	return 1;
277 
    }
278 
    return 0;
279 
}
280 
 
281 
static SEGP ctr_newseg(double x0, double y0, double x1, double y1, SEGP prev)
282 
{
283 
    SEGP seg = (SEGP)R_alloc(1, sizeof(SEG));
284 
    seg->x0 = x0;
285 
    seg->y0 = y0;
286 
    seg->x1 = x1;
287 
    seg->y1 = y1;
288 
    seg->next = prev;
289 
    return seg;
290 
}
291 
 
292 
static void ctr_swapseg(SEGP seg)
293 
{
294 
    double x, y;
295 
    x = seg->x0;
296 
    y = seg->y0;
297 
    seg->x0 = seg->x1;
298 
    seg->y0 = seg->y1;
299 
    seg->x1 = x;
300 
    seg->y1 = y;
301 
}
302 
 
13792 pd 303 
	/* ctr_segdir(): Determine the entry direction to the next cell */
1284 ihaka 304 
	/* and update the cell indices */
305 
 
13792 pd 306 
#define XMATCH(x0,x1) (fabs(x0-x1) == 0)
307 
#define YMATCH(y0,y1) (fabs(y0-y1) == 0)
1284 ihaka 308 
 
1839 ihaka 309 
static int ctr_segdir(double xend, double yend, double *x, double *y,
310 
		      int *i, int *j, int nx, int ny)
1284 ihaka 311 
{
1839 ihaka 312 
    if (YMATCH(yend, y[*j])) {
313 
	if (*j == 0)
314 
	    return 0;
1284 ihaka 315 
	*j = *j - 1;
316 
	return 3;
317 
    }
1839 ihaka 318 
    if (XMATCH(xend, x[*i])) {
319 
	if (*i == 0)
320 
	    return 0;
1284 ihaka 321 
	*i = *i - 1;
322 
	return 4;
323 
    }
1839 ihaka 324 
    if (YMATCH(yend, y[*j + 1])) {
325 
	if (*j >= ny - 1)
326 
	    return 0;
1284 ihaka 327 
	*j = *j + 1;
328 
	return 1;
329 
    }
1839 ihaka 330 
    if (XMATCH(xend, x[*i + 1])) {
331 
	if (*i >= nx - 1)
332 
	    return 0;
1284 ihaka 333 
	*i = *i + 1;
334 
	return 2;
335 
    }
336 
    return 0;
337 
}
338 
 
339 
/* Search seglist for a segment with endpoint (xend, yend). */
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 */
342 
/* is pointed to by seg and the updated segment list (with */
13792 pd 343 
/* the matched segment stripped) is returned by the funtion. */
1284 ihaka 344 
 
13792 pd 345 
static SEGP ctr_segupdate(double xend, double yend, int dir, Rboolean tail,
1839 ihaka 346 
			  SEGP seglist, SEGP* seg)
1284 ihaka 347 
{
1839 ihaka 348 
    if (seglist == NULL) {
1284 ihaka 349 
	*seg = NULL;
350 
	return NULL;
351 
    }
1839 ihaka 352 
    switch (dir) {
1284 ihaka 353 
    case 1:
354 
    case 3:
1839 ihaka 355 
	if (YMATCH(yend,seglist->y0)) {
356 
	    if (!tail)
357 
		ctr_swapseg(seglist);
1284 ihaka 358 
	    *seg = seglist;
359 
	    return seglist->next;
360 
	}
1839 ihaka 361 
	if (YMATCH(yend,seglist->y1)) {
362 
	    if (tail)
363 
		ctr_swapseg(seglist);
1284 ihaka 364 
	    *seg = seglist;
365 
	    return seglist->next;
366 
	}
367 
	break;
368 
    case 2:
369 
    case 4:
1839 ihaka 370 
	if (XMATCH(xend,seglist->x0)) {
371 
	    if (!tail)
372 
		ctr_swapseg(seglist);
1284 ihaka 373 
	    *seg = seglist;
374 
	    return seglist->next;
375 
	}
1839 ihaka 376 
	if (XMATCH(xend,seglist->x1)) {
377 
	    if (tail)
378 
		ctr_swapseg(seglist);
1284 ihaka 379 
	    *seg = seglist;
380 
	    return seglist->next;
381 
	}
382 
	break;
383 
    }
384 
    seglist->next = ctr_segupdate(xend, yend, dir, tail, seglist->next, seg);
385 
    return seglist;
386 
}
387 
 
7813 murrell 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().
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
9279 maechler 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
7813 murrell 394 
   different memory allocation mechanism.
395 
*/
396 
static SEXP labelList;
397 
 
7824 ripley 398 
static
13792 pd 399 
double distFromEdge(double *xxx, double *yyy, int iii, DevDesc *dd) {
17179 murrell 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]));
7813 murrell 402 
}
403 
 
7824 ripley 404 
static
13792 pd 405 
Rboolean useStart(double *xxx, double *yyy, int ns, DevDesc *dd) {
7813 murrell 406 
    if (distFromEdge(xxx, yyy, 0, dd) < distFromEdge(xxx, yyy, ns-1, dd))
13792 pd 407 
	return TRUE;
9279 maechler 408 
    else
13792 pd 409 
	return FALSE;
7813 murrell 410 
}
411 
 
7824 ripley 412 
static
9279 maechler 413 
int findGapUp(double *xxx, double *yyy, int ns, double labelDistance,
7813 murrell 414 
	      DevDesc *dd) {
415 
    double dX, dY;
416 
    double dXC, dYC;
417 
    double distanceSum = 0;
418 
    int n = 0;
419 
    int jjj = 1;
420 
    while ((jjj < ns) && (distanceSum < labelDistance)) {
421 
	/* Find a gap big enough for the label
422 
	   use several segments if necessary
423 
	*/
424 
	dX = xxx[jjj] - xxx[jjj - n - 1];
425 
	dY = yyy[jjj] - yyy[jjj - n - 1];
426 
	dXC = GConvertXUnits(dX, USER, INCHES, dd);
427 
	dYC = GConvertYUnits(dY, USER, INCHES, dd);
9280 maechler 428 
	distanceSum = hypot(dXC, dYC);
7813 murrell 429 
	jjj = (jjj + 1);
430 
	n += 1;
431 
    }
432 
    if (distanceSum < labelDistance)
433 
	return 0;
434 
    else
435 
	return n;
436 
}
437 
 
7824 ripley 438 
static
7813 murrell 439 
int findGapDown(double *xxx, double *yyy, int ns, double labelDistance,
440 
		DevDesc *dd) {
441 
    double dX, dY;
442 
    double dXC, dYC;
443 
    double distanceSum = 0;
444 
    int n = 0;
445 
    int jjj = ns - 2;
446 
    while ((jjj > -1) && (distanceSum < labelDistance)) {
447 
	/* Find a gap big enough for the label
448 
	   use several segments if necessary
449 
	*/
450 
	dX = xxx[jjj] - xxx[jjj + n + 1];
451 
	dY = yyy[jjj] - yyy[jjj + n + 1];
452 
	dXC = GConvertXUnits(dX, USER, INCHES, dd);
453 
	dYC = GConvertYUnits(dY, USER, INCHES, dd);
9280 maechler 454 
	distanceSum = hypot(dXC, dYC);
9279 maechler 455 
	jjj--;
7813 murrell 456 
	n -= 1;
457 
    }
458 
    if (distanceSum < labelDistance)
459 
	return 0;
460 
    else
461 
	return -n;
462 
}
9279 maechler 463 
 
21061 murrell 464 
/*
465 
 * Generate a list of segments for a single level
466 
 */
467 
static SEGP* contourLines(double *x, int nx, double *y, int ny,
468 
			 double *z, double zc, double atom)
469 
{
470 
    double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];
471 
    int i, j, k, l, m, nacode;
472 
    SEGP seglist;
473 
    SEGP *segmentDB;
474 
    /* Initialize the segment data base */
475 
    /* Note we must be careful about resetting */
476 
    /* the top of the stack, otherwise we run out of */
477 
    /* memory after a sequence of displaylist replays */
478 
    /*
479 
     * This reset is done out in GEcontourLines
480 
     */
481 
    segmentDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));
482 
    for (i = 0; i < nx; i++)
483 
	for (j = 0; j < ny; j++)
484 
	    segmentDB[i + j * nx] = NULL;
485 
    for (i = 0; i < nx - 1; i++) {
486 
	xl = x[i];
487 
	xh = x[i + 1];
488 
	for (j = 0; j < ny - 1; j++) {
489 
	    yl = y[j];
490 
	    yh = y[j + 1];
491 
	    k = i + j * nx;
492 
	    zll = z[k];
493 
	    zhl = z[k + 1];
494 
	    zlh = z[k + nx];
495 
	    zhh = z[k + nx + 1];
496 
 
497 
	    /* If the value at a corner is exactly equal to a contour level,
498 
	     * change that value by a tiny amount */
499 
 
500 
	    if (zll == zc) zll += atom;
501 
	    if (zhl == zc) zhl += atom;
502 
	    if (zlh == zc) zlh += atom;
503 
	    if (zhh == zc) zhh += atom;
504 
#ifdef DEBUG_contour
505 
	    /* Haven't seen this happening (MM): */
506 
	    if (zll == zc) REprintf(" [%d,%d] ll: %g\n",i,j, zll);
507 
	    if (zhl == zc) REprintf(" [%d,%d] hl: %g\n",i,j, zhl);
508 
	    if (zlh == zc) REprintf(" [%d,%d] lh: %g\n",i,j, zlh);
509 
	    if (zhh == zc) REprintf(" [%d,%d] hh: %g\n",i,j, zhh);
510 
#endif
511 
	    /* Check for intersections with sides */
512 
 
513 
	    nacode = 0;
514 
	    if (R_FINITE(zll)) nacode += 1;
515 
	    if (R_FINITE(zhl)) nacode += 2;
516 
	    if (R_FINITE(zlh)) nacode += 4;
517 
	    if (R_FINITE(zhh)) nacode += 8;
518 
 
519 
	    k = 0;
520 
	    switch (nacode) {
521 
	    case 15:
522 
		if (ctr_intersect(zll, zhl, zc, &f)) {
523 
		    xx[k] = xl + f * (xh - xl);
524 
		    yy[k] = yl; k++;
525 
		}
526 
		if (ctr_intersect(zll, zlh, zc, &f)) {
527 
		    yy[k] = yl + f * (yh - yl);
528 
		    xx[k] = xl; k++;
529 
		}
530 
		if (ctr_intersect(zhl, zhh, zc, &f)) {
531 
		    yy[k] = yl + f * (yh - yl);
532 
		    xx[k] = xh; k++;
533 
		}
534 
		if (ctr_intersect(zlh, zhh, zc, &f)) {
535 
		    xx[k] = xl + f * (xh - xl);
536 
		    yy[k] = yh; k++;
537 
		}
538 
		break;
539 
	    case 14:
540 
		if (ctr_intersect(zhl, zhh, zc, &f)) {
541 
		    yy[k] = yl + f * (yh - yl);
542 
		    xx[k] = xh; k++;
543 
		}
544 
		if (ctr_intersect(zlh, zhh, zc, &f)) {
545 
		    xx[k] = xl + f * (xh - xl);
546 
		    yy[k] = yh; k++;
547 
		}
548 
		if (ctr_intersect(zlh, zhl, zc, &f)) {
549 
		    xx[k] = xl + f * (xh - xl);
550 
		    yy[k] = yh + f * (yl - yh);
551 
		    k++;
552 
		}
553 
		break;
554 
	    case 13:
555 
		if (ctr_intersect(zll, zlh, zc, &f)) {
556 
		    yy[k] = yl + f * (yh - yl);
557 
		    xx[k] = xl; k++;
558 
		}
559 
		if (ctr_intersect(zlh, zhh, zc, &f)) {
560 
		    xx[k] = xl + f * (xh - xl);
561 
		    yy[k] = yh; k++;
562 
		}
563 
		if (ctr_intersect(zll, zhh, zc, &f)) {
564 
		    xx[k] = xl + f * (xh - xl);
565 
		    yy[k] = yl + f * (yh - yl);
566 
		    k++;
567 
		}
568 
		break;
569 
	    case 11:
570 
		if (ctr_intersect(zhl, zhh, zc, &f)) {
571 
		    yy[k] = yl + f * (yh - yl);
572 
		    xx[k] = xh; k++;
573 
		}
574 
		if (ctr_intersect(zll, zhl, zc, &f)) {
575 
		    xx[k] = xl + f * (xh - xl);
576 
		    yy[k] = yl; k++;
577 
		}
578 
		if (ctr_intersect(zll, zhh, zc, &f)) {
579 
		    xx[k] = xl + f * (xh - xl);
580 
		    yy[k] = yl + f * (yh - yl);
581 
		    k++;
582 
		}
583 
		break;
584 
	    case 7:
585 
		if (ctr_intersect(zll, zlh, zc, &f)) {
586 
		    yy[k] = yl + f * (yh - yl);
587 
		    xx[k] = xl; k++;
588 
		}
589 
		if (ctr_intersect(zll, zhl, zc, &f)) {
590 
		    xx[k] = xl + f * (xh - xl);
591 
		    yy[k] = yl; k++;
592 
		}
593 
		if (ctr_intersect(zlh, zhl, zc, &f)) {
594 
		    xx[k] = xl + f * (xh - xl);
595 
		    yy[k] = yh + f * (yl - yh);
596 
		    k++;
597 
		}
598 
		break;
599 
	    }
600 
 
601 
	    /* We now have k(=2,4) endpoints */
602 
	    /* Decide which to join */
603 
 
604 
	    seglist = NULL;
605 
 
606 
	    if (k > 0) {
607 
		if (k == 2) {
608 
		    seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
609 
		}
610 
		else if (k == 4) {
611 
		    for (k = 3; k >= 1; k--) {
612 
			m = k;
613 
			xl = xx[k];
614 
			for (l = 0; l < k; l++) {
615 
			    if (xx[l] > xl) {
616 
				xl = xx[l];
617 
				m = l;
618 
			    }
619 
			}
620 
			if (m != k) {
621 
			    xl = xx[k];
622 
			    yl = yy[k];
623 
			    xx[k] = xx[m];
624 
			    yy[k] = yy[m];
625 
			    xx[m] = xl;
626 
			    yy[m] = yl;
627 
			}
628 
		    }
629 
		    seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
630 
		    seglist = ctr_newseg(xx[2], yy[2], xx[3], yy[3], seglist);
631 
		}
632 
		else error("k != 2 or 4");
633 
	    }
634 
	    segmentDB[i + j * nx] = seglist;
635 
	}
636 
    }
637 
    return segmentDB;
638 
}
639 
 
640 
/* maximal number of line segments of one contour segment:
641 
 * for preventing infinite loops -- shouldn't be needed --> warning */
642 
#define MAX_ns 25000
643 
 
644 
#define CONTOUR_LIST_STEP 100
645 
#define CONTOUR_LIST_LEVEL 0
646 
#define CONTOUR_LIST_X 1
647 
#define CONTOUR_LIST_Y 2
648 
 
649 
static SEXP growList(SEXP oldlist) {
650 
    int i, len;
651 
    SEXP templist;
652 
    len = LENGTH(oldlist);
653 
    templist = PROTECT(allocVector(VECSXP, len + CONTOUR_LIST_STEP));
654 
    for (i=0; i<len; i++)
655 
	SET_VECTOR_ELT(templist, i, VECTOR_ELT(oldlist, i));
656 
    UNPROTECT(1);
657 
    return templist;
658 
}
659 
 
660 
/*
661 
 * Store the list of segments for a single level in the SEXP
662 
 * list that will be returned to the user
663 
 */
664 
int addContourLines(double *x, int nx, double *y, int ny,
665 
		     double *z, double zc, double atom,
666 
		     SEGP* segmentDB, int nlines, SEXP container)
667 
{
668 
    double xend, yend;
669 
    int i, ii, j, jj, ns, ns2, dir, nc;
670 
    SEGP seglist, seg, s, start, end;
21594 murrell 671 
    SEXP ctr, level, xsxp, ysxp, names;
21061 murrell 672 
    /* Begin following contours. */
673 
    /* 1. Grab a segment */
674 
    /* 2. Follow its tail */
675 
    /* 3. Follow its head */
676 
    /* 4. Save the contour */
677 
    for (i = 0; i < nx - 1; i++)
678 
	for (j = 0; j < ny - 1; j++) {
679 
	    while ((seglist = segmentDB[i + j * nx])) {
680 
		ii = i; jj = j;
681 
		start = end = seglist;
682 
		segmentDB[i + j * nx] = seglist->next;
683 
		xend = seglist->x1;
684 
		yend = seglist->y1;
685 
		while ((dir = ctr_segdir(xend, yend, x, y,
686 
					 &ii, &jj, nx, ny))) {
687 
		    segmentDB[ii + jj * nx]
688 
			= ctr_segupdate(xend, yend, dir, TRUE,/* = tail */
689 
					segmentDB[ii + jj * nx], &seg);
690 
		    if (!seg) break;
691 
		    end->next = seg;
692 
		    end = seg;
693 
		    xend = end->x1;
694 
		    yend = end->y1;
695 
		}
696 
		end->next = NULL; /* <<< new for 1.2.3 */
697 
		ii = i; jj = j;
698 
		xend = seglist->x0;
699 
		yend = seglist->y0;
700 
		while ((dir = ctr_segdir(xend, yend, x, y,
701 
					 &ii, &jj, nx, ny))) {
702 
		    segmentDB[ii + jj * nx]
703 
			= ctr_segupdate(xend, yend, dir, FALSE,/* ie. head */
704 
					segmentDB[ii+jj*nx], &seg);
705 
		    if (!seg) break;
706 
		    seg->next = start;
707 
		    start = seg;
708 
		    xend = start->x0;
709 
		    yend = start->y0;
710 
		}
711 
 
712 
		/* ns := #{segments of polyline} -- need to allocate */
713 
		s = start;
714 
		ns = 0;
715 
		/* MAX_ns: prevent inf.loop (shouldn't be needed) */
716 
		while (s && ns < MAX_ns) {
717 
		    ns++;
718 
		    s = s->next;
719 
		}
720 
		if(ns == MAX_ns)
721 
		    warning("contour(): circular/long seglist -- bug.report()!");
722 
 
723 
		/* countour midpoint : use for labelling sometime (not yet!) */
724 
		if (ns > 3) ns2 = ns/2; else ns2 = -1;
725 
 
726 
		/*
727 
		 * "write" the contour locations into the list of contours
728 
		 */
729 
		ctr = PROTECT(allocVector(VECSXP, 3));
730 
		level = PROTECT(allocVector(REALSXP, 1));
731 
		xsxp = PROTECT(allocVector(REALSXP, ns + 1));
732 
		ysxp = PROTECT(allocVector(REALSXP, ns + 1));
733 
		REAL(level)[0] = zc;
734 
		SET_VECTOR_ELT(ctr, CONTOUR_LIST_LEVEL, level);
735 
		s = start;
736 
		REAL(xsxp)[0] = s->x0;
737 
		REAL(ysxp)[0] = s->y0;
738 
		ns = 1;
739 
		while (s->next && ns < MAX_ns) {
740 
		    s = s->next;
741 
		    REAL(xsxp)[ns] = s->x0;
742 
		    REAL(ysxp)[ns++] = s->y0;
743 
		}
744 
		REAL(xsxp)[ns] = s->x1;
745 
		REAL(ysxp)[ns] = s->y1;
746 
		SET_VECTOR_ELT(ctr, CONTOUR_LIST_X, xsxp);
747 
		SET_VECTOR_ELT(ctr, CONTOUR_LIST_Y, ysxp);
748 
		/*
21594 murrell 749 
		 * Set the names attribute for the contour
750 
		 * So that users can extract components using
751 
		 * meaningful names
752 
		 */
753 
		PROTECT(names = allocVector(STRSXP, 3));
754 
		SET_STRING_ELT(names, 0, mkChar("level"));
755 
		SET_STRING_ELT(names, 1, mkChar("x"));
756 
		SET_STRING_ELT(names, 2, mkChar("y"));
757 
		setAttrib(ctr, R_NamesSymbol, names);
758 
		/*
21061 murrell 759 
		 * We're about to add another line to the list ...
760 
		 */
761 
		nlines += 1;
762 
		nc = LENGTH(VECTOR_ELT(container, 0));
763 
		if (nlines == nc)
764 
		    /* Where does this get UNPROTECTed? */
765 
		    SET_VECTOR_ELT(container, 0,
766 
				   growList(VECTOR_ELT(container, 0)));
767 
		SET_VECTOR_ELT(VECTOR_ELT(container, 0), nlines - 1, ctr);
21594 murrell 768 
		UNPROTECT(5);
21061 murrell 769 
	    }
770 
	}
771 
    return nlines;
772 
}
773 
 
774 
/*
775 
 * Given nx x values, ny y values, nx*ny z values,
776 
 * and nl cut-values in z ...
777 
 * ... produce a list of contour lines:
778 
 *   list of sub-lists
779 
 *     sub-list = x vector, y vector, and cut-value.
780 
 */
781 
SEXP GEcontourLines(double *x, int nx, double *y, int ny,
782 
		    double *z, double *levels, int nl,
783 
		    GEDevDesc *dd)
784 
{
785 
    char *vmax;
786 
    int i, nlines, len;
787 
    double atom, zmin, zmax;
788 
    SEGP* segmentDB;
789 
    SEXP container, mainlist, templist;
790 
    /*
791 
     * "tie-breaker" values
792 
     */
793 
    zmin = DBL_MAX;
794 
    zmax = DBL_MIN;
795 
    for (i = 0; i < nx * ny; i++)
796 
	if (R_FINITE(z[i])) {
797 
	    if (zmax < z[i]) zmax =  z[i];
798 
	    if (zmin > z[i]) zmin =  z[i];
799 
	}
800 
 
801 
    if (zmin >= zmax) {
802 
	if (zmin == zmax)
803 
	    warning("all z values are equal");
804 
	else
805 
	    warning("all z values are NA");
806 
	return R_NilValue;
807 
    }
808 
    /* change to 1e-3, reconsidered because of PR#897
809 
     * but 1e-7, and even  2*DBL_EPSILON do not prevent inf.loop in contour().
810 
     * maybe something like   16 * DBL_EPSILON * (..).
811 
     * see also MAX_ns above */
812 
    atom = 1e-3 * (zmax - zmin);
813 
    /*
814 
     * Create a "container" which is a list with only 1 element.
815 
     * The element is the list of lines that will be built up.
816 
     * I create the container because this allows me to PROTECT
817 
     * the container once here and then UNPROTECT it at the end of
818 
     * this function and, as long as I always work with
819 
     * VECTOR_ELT(container, 0) and SET_VECTOR_ELT(container, 0)
820 
     * in functions called from here, I don't need to worry about
821 
     * protectin the list that I am building up.
822 
     * Why bother?  Because the list I am building can potentially
823 
     * grow and it's awkward to get the PROTECTs/UNPROTECTs right
824 
     * when you're in a loop and growing a list.
825 
     */
826 
    container = PROTECT(allocVector(VECSXP, 1));
827 
    /*
828 
     * Create "large" list (will trim excess at the end if necesary)
829 
     */
830 
    SET_VECTOR_ELT(container, 0, allocVector(VECSXP, CONTOUR_LIST_STEP));
831 
    nlines = 0;
832 
    /*
833 
     * Add lines for each contour level
834 
     */
835 
    for (i = 0; i < nl; i++) {
836 
	/*
837 
	 * The vmaxget/set is to manage the memory that gets
838 
	 * R_alloc'ed in the creation of the segmentDB structure
839 
	 */
840 
	vmax = vmaxget();
841 
	/*
842 
	 * Generate a segment database
843 
	 */
844 
	segmentDB = contourLines(x, nx, y, ny, z, levels[i], atom);
845 
	/*
846 
	 * Add lines to the list based on the segment database
847 
	 */
848 
	nlines = addContourLines(x, nx, y, ny, z, levels[i],
849 
				 atom, segmentDB, nlines,
850 
				 container);
851 
	vmaxset(vmax);
852 
    }
853 
    /*
854 
     * Trim the list of lines to the appropriate length.
855 
     */
856 
    len = LENGTH(VECTOR_ELT(container, 0));
857 
    if (nlines < len) {
858 
	mainlist = VECTOR_ELT(container, 0);
859 
	templist = PROTECT(allocVector(VECSXP, nlines));
860 
	for (i=0; i<nlines; i++)
861 
	    SET_VECTOR_ELT(templist, i, VECTOR_ELT(mainlist, i));
862 
	mainlist = templist;
863 
	UNPROTECT(1);  /* UNPROTECT templist */
864 
    } else
865 
	mainlist = VECTOR_ELT(container, 0);
866 
    UNPROTECT(1);  /* UNPROTECT container */
867 
    return mainlist;
868 
}
869 
 
870 
SEXP GEdrawContourLines();
871 
 
27215 murrell 872 
SEXP do_contourLines(SEXP call, SEXP op, SEXP args, SEXP env)
873 
{
874 
    SEXP oargs, c, x, y, z;
875 
    int nx, ny, nc;
876 
    GEDevDesc *dd = GEcurrentDevice();
877 
 
878 
    oargs = args;
879 
 
880 
    x = CAR(args);
881 
    internalTypeCheck(call, x, REALSXP);
882 
    nx = LENGTH(x);
883 
    args = CDR(args);
884 
 
885 
    y = CAR(args);
886 
    internalTypeCheck(call, y, REALSXP);
887 
    ny = LENGTH(y);
888 
    args = CDR(args);
889 
 
890 
    z = CAR(args);
891 
    internalTypeCheck(call, z, REALSXP);
892 
    args = CDR(args);
893 
 
894 
    /* levels */
895 
    c = CAR(args);
896 
    internalTypeCheck(call, c, REALSXP);
897 
    nc = LENGTH(c);
898 
    args = CDR(args);
899 
 
900 
    return GEcontourLines(REAL(x), nx, REAL(y), ny, REAL(z), REAL(c), nc, dd);
901 
}
902 
 
903 
/*
904 
 * The *base graphics* function contour() and the *general base*
905 
 * function contourLines() use the same code to generate contour lines
906 
 * (i.e., the function contourLines())
907 
 *
908 
 * I had a look at extracting the code that draws the labels
909 
 * into a *general base* function
910 
 * (e.g., into some sort of labelLines() function),
911 
 * but the code is too base-graphics-specific (e.g., one of the
912 
 * labelling methods seeks the location closest to the edge of the
913 
 * plotting region) so I've left it alone for now.
914 
 *
915 
 * This does mean that the contourLines() function is part of the
916 
 * graphics engine, but the contour() function is part of the
917 
 * base graphics system.
918 
 */
919 
 
13792 pd 920 
static void contour(SEXP x, int nx, SEXP y, int ny, SEXP z,
921 
		    double zc,
8213 murrell 922 
		    SEXP labels, int cnum,
10886 maechler 923 
		    Rboolean drawLabels, int method,
924 
		    Rboolean vectorFonts, int typeface, int fontindex,
925 
		    double atom, DevDesc *dd)
1284 ihaka 926 
{
13792 pd 927 
/* draw a contour for one given contour level `zc' */
928 
 
15321 luke 929 
    char *vmax;
930 
 
1284 ihaka 931 
    double xend, yend;
27215 murrell 932 
    int i, ii, j, jj, ns, ns2, dir;
1284 ihaka 933 
    SEGP seglist, seg, s, start, end;
934 
    double *xxx, *yyy;
935 
 
7813 murrell 936 
    double variance, dX, dY, deltaX, deltaY;
937 
    double dXC, dYC, deltaXC, deltaYC;
11067 maechler 938 
    int range=0, indx=0, n; /* -Wall */
7813 murrell 939 
    double lowestVariance;
940 
    double squareSum, sum;
941 
    int iii, jjj;
942 
    double distanceSum, labelDistance, avgGradient;
943 
    int zeroCount;
944 
    char buffer[255];
7824 ripley 945 
    double avg;
7813 murrell 946 
    int result;
947 
    double ux, uy, vx, vy;
7824 ripley 948 
    double xStart, yStart;
7813 murrell 949 
    double dx, dy, dxy;
950 
    double labelHeight;
8872 pd 951 
    SEXP label1 = PROTECT(allocVector(REALSXP, 8));
7813 murrell 952 
    SEXP label2;
8213 murrell 953 
    SEXP lab;
10886 maechler 954 
    Rboolean gotLabel = FALSE;
955 
    Rboolean ddl;/* Don't draw label -- currently unused, i.e. always FALSE*/
7813 murrell 956 
 
13792 pd 957 
#ifdef DEBUG_contour
958 
    Rprintf("contour(lev = %g):\n", zc);
959 
#endif
1284 ihaka 960 
 
27215 murrell 961 
    vmax = vmaxget();
962 
    ctr_SegDB = contourLines(REAL(x), nx, REAL(y), ny, REAL(z), zc, atom);
27811 ripley 963 
    /* we need to keep ctr_SegDB available, so vmaxset(vmax); was wrong */
1284 ihaka 964 
 
965 
    /* The segment database is now assembled. */
966 
    /* Begin following contours. */
967 
    /* 1. Grab a segment */
968 
    /* 2. Follow its tail */
969 
    /* 3. Follow its head */
970 
    /* 4. Draw the contour */
971 
 
1839 ihaka 972 
    for (i = 0; i < nx - 1; i++)
10886 maechler 973 
      for (j = 0; j < ny - 1; j++) {
974 
	while ((seglist = ctr_SegDB[i + j * nx])) {
975 
	    ii = i; jj = j;
976 
	    start = end = seglist;
977 
	    ctr_SegDB[i + j * nx] = seglist->next;
978 
	    xend = seglist->x1;
979 
	    yend = seglist->y1;
980 
	    while ((dir = ctr_segdir(xend, yend, REAL(x), REAL(y),
981 
				     &ii, &jj, nx, ny))) {
982 
		ctr_SegDB[ii + jj * nx]
13792 pd 983 
		    = ctr_segupdate(xend, yend, dir, TRUE,/* = tail */
10886 maechler 984 
				    ctr_SegDB[ii + jj * nx], &seg);
985 
		if (!seg) break;
986 
		end->next = seg;
987 
		end = seg;
988 
		xend = end->x1;
989 
		yend = end->y1;
990 
	    }
13792 pd 991 
	    end->next = NULL; /* <<< new for 1.2.3 */
10886 maechler 992 
	    ii = i; jj = j;
993 
	    xend = seglist->x0;
994 
	    yend = seglist->y0;
995 
	    while ((dir = ctr_segdir(xend, yend, REAL(x), REAL(y),
996 
				     &ii, &jj, nx, ny))) {
997 
		ctr_SegDB[ii + jj * nx]
13792 pd 998 
		    = ctr_segupdate(xend, yend, dir, FALSE,/* ie. head */
10886 maechler 999 
				    ctr_SegDB[ii+jj*nx], &seg);
1000 
		if (!seg) break;
1001 
		seg->next = start;
1002 
		start = seg;
1003 
		xend = start->x0;
1004 
		yend = start->y0;
1005 
	    }
13792 pd 1006 
 
1007 
	    /* ns := #{segments of polyline} -- need to allocate */
10886 maechler 1008 
	    s = start;
1009 
	    ns = 0;
13792 pd 1010 
	    /* MAX_ns: prevent inf.loop (shouldn't be needed) */
1011 
	    while (s && ns < MAX_ns) {
10886 maechler 1012 
		ns++;
1013 
		s = s->next;
1014 
	    }
13792 pd 1015 
	    if(ns == MAX_ns)
1016 
		warning("contour(): circular/long seglist -- bug.report()!");
1284 ihaka 1017 
 
13792 pd 1018 
	    /* countour midpoint : use for labelling sometime (not yet!) */
1019 
	    if (ns > 3) ns2 = ns/2; else ns2 = -1;
1284 ihaka 1020 
 
15321 luke 1021 
	    vmax = vmaxget();
1022 
	    xxx = (double *) R_alloc(ns + 1, sizeof(double));
1023 
	    yyy = (double *) R_alloc(ns + 1, sizeof(double));
13792 pd 1024 
	    /* now have the space, go through again: */
1025 
	    s = start;
10886 maechler 1026 
	    ns = 0;
1027 
	    xxx[ns] = s->x0;
1028 
	    yyy[ns++] = s->y0;
13792 pd 1029 
	    while (s->next && ns < MAX_ns) {
10886 maechler 1030 
		s = s->next;
1284 ihaka 1031 
		xxx[ns] = s->x0;
1032 
		yyy[ns++] = s->y0;
10886 maechler 1033 
	    }
1034 
	    xxx[ns] = s->x1;
1035 
	    yyy[ns++] = s->y1;
13792 pd 1036 
#ifdef DEBUG_contour
1037 
	    Rprintf("  [%2d,%2d]: (x,y)[1:%d] = ", i,j, ns);
1038 
	    if(ns >= 5)
1039 
		Rprintf(" (%g,%g), (%g,%g), ..., (%g,%g)\n",
1040 
			xxx[0],yyy[0], xxx[1],yyy[1], xxx[ns-1],yyy[ns-1]);
1041 
	    else
1042 
		for(iii = 0; iii < ns; iii++)
1043 
		    Rprintf(" (%g,%g)%s", xxx[iii],yyy[iii],
1044 
			    (iii < ns-1) ? "," : "\n");
1045 
#endif
1046 
 
10886 maechler 1047 
	    GMode(1, dd);
7813 murrell 1048 
 
13792 pd 1049 
	    if (drawLabels) {
1050 
		/* If user supplied labels, use i'th one of them
1051 
		   Otherwise stringify the z-value of the contour */
1052 
		buffer[0] = ' ';
1053 
		if (!isNull(labels)) {
1054 
		    int numl = length(labels);
1055 
		    strcpy(&buffer[1], CHAR(STRING_ELT(labels, cnum % numl)));
1056 
		}
1057 
		else {
1058 
		    PROTECT(lab = allocVector(REALSXP, 1));
1059 
		    REAL(lab)[0] = zc;
1060 
		    lab = labelformat(lab);
1061 
		    strcpy(&buffer[1], CHAR(STRING_ELT(lab, 0)));
1062 
		    UNPROTECT(1);
1063 
		}
1064 
		buffer[strlen(buffer)+1] = '\0';
1065 
		buffer[strlen(buffer)] = ' ';
7813 murrell 1066 
 
13792 pd 1067 
		if (vectorFonts) {
1068 
		    /* 1, 1 => sans serif, basic font */
1069 
		    labelDistance =
1070 
			GVStrWidth((unsigned char *)buffer, typeface, fontindex,
1071 
				   INCHES, dd);
1072 
		    labelHeight =
1073 
			GVStrHeight((unsigned char *)buffer, typeface, fontindex,
1074 
				    INCHES, dd);
1075 
		}
1076 
		else {
1077 
		    labelDistance = GStrWidth(buffer, INCHES, dd);
1078 
		    labelHeight = GStrHeight(buffer, INCHES, dd);
1079 
		}
7813 murrell 1080 
 
13792 pd 1081 
		if (labelDistance > 0) {
1082 
		    /* Try to find somewhere to draw the label */
1083 
		    switch (method) {
1084 
		    case 0: /* draw label at one end of contour
1085 
			       overwriting contour line
1086 
			    */
1087 
			if (useStart(xxx, yyy, ns, dd) )
1088 
			    indx = 0;
1089 
			else
1090 
			    indx = ns - 1;
1091 
			break;
1092 
		    case 1: /* draw label at one end of contour
1093 
			       embedded in contour
1094 
			       no overlapping labels
1095 
			    */
11067 maechler 1096 
			indx = 0;
13792 pd 1097 
			range = 0;
1098 
			gotLabel = FALSE;
1099 
			if (useStart(xxx, yyy, ns, dd)) {
1100 
			    iii = 0;
1101 
			    n = findGapUp(xxx, yyy, ns, labelDistance, dd);
1102 
			}
1103 
			else {
1104 
			    n = findGapDown(xxx, yyy, ns, labelDistance, dd);
1105 
			    iii = ns - n - 1;
1106 
			}
1107 
			if (n > 0) {
1108 
			    /** Find 4 corners of label extents **/
1109 
			    FindCorners(labelDistance, labelHeight, label1,
1110 
					xxx[iii], yyy[iii],
1111 
					xxx[iii+n], yyy[iii+n], dd);
10886 maechler 1112 
 
13792 pd 1113 
			    /** Test corners for intersection with previous labels **/
1114 
			    label2 = labelList;
1115 
			    result = 0;
1116 
			    while ((result == 0) && (label2 != R_NilValue)) {
1117 
				result = TestLabelIntersection(label1, CAR(label2));
1118 
				label2 = CDR(label2);
1119 
			    }
7813 murrell 1120 
			    if (result == 0) {
13792 pd 1121 
				result = LabelInsideWindow(label1, dd);
1122 
				if (result == 0) {
1123 
				    indx = iii;
1124 
				    range = n;
1125 
				    gotLabel = TRUE;
1126 
				}
7813 murrell 1127 
			    }
1128 
			}
13792 pd 1129 
			break;
1130 
		    case 2: /* draw label on flattest portion of contour
1131 
			       embedded in contour line
1132 
			       no overlapping labels
1133 
			    */
1134 
			/* Look for flatest sequence of contour gradients */
1135 
			lowestVariance = 9999999;   /* A large number */
1136 
			indx = 0;
1137 
			range = 0;
1138 
			gotLabel = FALSE;
1139 
			for (iii = 0; iii < ns; iii++) {
1140 
			    distanceSum = 0;
1141 
			    avgGradient = 0;
1142 
			    squareSum = 0;
1143 
			    sum = 0;
1144 
			    n = 0;
1145 
			    zeroCount = 0;
1146 
			    jjj = (iii + 1);
1147 
			    while ((jjj < ns-1) &&
1148 
				   (distanceSum < labelDistance)) {
9279 maechler 1149 
 
7813 murrell 1150 
				/* Find a gap big enough for the label
1151 
				   use several segments if necessary
1152 
				*/
13792 pd 1153 
				dX = xxx[jjj] - xxx[jjj - n - 1];
1154 
				dY = yyy[jjj] - yyy[jjj - n - 1];
1155 
				dXC = GConvertXUnits(dX, USER, INCHES, dd);
1156 
				dYC = GConvertYUnits(dY, USER, INCHES, dd);
1157 
				distanceSum = hypot(dXC, dYC);
9279 maechler 1158 
 
7813 murrell 1159 
				/* Calculate the variance of the gradients
1160 
				   of the segments that will make way for the
1161 
				   label
1162 
				*/
13792 pd 1163 
				deltaX = xxx[jjj] - xxx[jjj - 1];
1164 
				deltaY = yyy[jjj] - yyy[jjj - 1];
1165 
				deltaXC = GConvertXUnits(deltaX, USER, INCHES, dd);
1166 
				deltaYC = GConvertYUnits(deltaY, USER, INCHES, dd);
1167 
				if (deltaX == 0) {deltaX = 1;}
1168 
				avgGradient += (deltaY/deltaX);
1169 
				squareSum += avgGradient * avgGradient;
1170 
				jjj = (jjj + 1);
1171 
				n += 1;
1172 
			    }
1173 
			    if (distanceSum < labelDistance)
1174 
				break;
9279 maechler 1175 
 
13792 pd 1176 
			    /** Find 4 corners of label extents **/
1177 
			    FindCorners(labelDistance, labelHeight, label1,
1178 
					xxx[iii], yyy[iii],
1179 
					xxx[iii+n], yyy[iii+n], dd);
9279 maechler 1180 
 
13792 pd 1181 
			    /** Test corners for intersection with previous labels **/
1182 
			    label2 = labelList;
1183 
			    result = 0;
1184 
			    while ((result == 0) && (label2 != R_NilValue)) {
1185 
				result = TestLabelIntersection(label1, CAR(label2));
1186 
				label2 = CDR(label2);
7813 murrell 1187 
			    }
13792 pd 1188 
			    if (result == 0)
1189 
				result = LabelInsideWindow(label1, dd);
1190 
			    if (result == 0) {
1191 
				variance = (squareSum - (avgGradient * avgGradient) / n) / n;
1192 
				avgGradient /= n;
1193 
				if (variance < lowestVariance) {
1194 
				    lowestVariance = variance;
1195 
				    indx = iii;
1196 
				    range = n;
1197 
				    avg = avgGradient;
1198 
				}
1199 
			    }
1200 
			    if (lowestVariance < 9999999)
1201 
				gotLabel = TRUE;
7813 murrell 1202 
			}
13792 pd 1203 
		    } /* switch (method) */
1204 
 
1205 
		    if (method == 0) {
1206 
			GPolyline(ns, xxx, yyy, USER, dd);
1207 
			if (vectorFonts)
1208 
			    GVText(xxx[indx], yyy[indx], USER, buffer,
1209 
				   typeface, fontindex,
1210 
				   .5, .5, 0, dd);
1211 
			else
1212 
			    GText(xxx[indx], yyy[indx], USER, buffer,
1213 
				  .5, .5, 0, dd);
10886 maechler 1214 
		    }
13792 pd 1215 
		    else {
1216 
			for (iii = 0; iii < indx; iii++)
1217 
			    GLine(xxx[iii], yyy[iii],
1218 
				  xxx[iii+1], yyy[iii+1], USER, dd);
1219 
			for (iii = indx+range; iii < ns - 1; iii++)
1220 
			    GLine(xxx[iii], yyy[iii],
1221 
				  xxx[iii+1], yyy[iii+1], USER, dd);
9279 maechler 1222 
 
13792 pd 1223 
			if (gotLabel) {
1224 
			    /* find which plot edge we are closest to */
1225 
			    int closest; /* 0 = indx,  1 = indx+range */
1226 
			    double dx1, dx2, dy1, dy2, dmin;
17179 murrell 1227 
			    dx1 = fmin2((xxx[indx] - Rf_gpptr(dd)->usr[0]),
1228 
					(Rf_gpptr(dd)->usr[1] - xxx[indx]));
1229 
			    dx2 = fmin2((Rf_gpptr(dd)->usr[1] - xxx[indx+range]),
1230 
					(xxx[indx+range] - Rf_gpptr(dd)->usr[0]));
13792 pd 1231 
			    if (dx1 < dx2) {
1232 
				closest = 0;
1233 
				dmin = dx1;
1234 
			    } else {
1235 
				closest = 1;
1236 
				dmin = dx2;
1237 
			    }
17179 murrell 1238 
			    dy1 = fmin2((yyy[indx] - Rf_gpptr(dd)->usr[2]),
1239 
					(Rf_gpptr(dd)->usr[3] - yyy[indx]));
13792 pd 1240 
			    if (closest && (dy1 < dmin)) {
1241 
				closest = 0;
1242 
				dmin = dy1;
1243 
			    } else if (dy1 < dmin)
1244 
				dmin = dy1;
17179 murrell 1245 
			    dy2 = fmin2((Rf_gpptr(dd)->usr[3] - yyy[indx+range]),
1246 
					(yyy[indx+range] - Rf_gpptr(dd)->usr[2]));
13792 pd 1247 
			    if (!closest && (dy2 < dmin))
1248 
				closest = 1;
7813 murrell 1249 
 
13792 pd 1250 
			    dx = GConvertXUnits(xxx[indx+range] - xxx[indx],
1251 
						USER, INCHES, dd);
1252 
			    dy = GConvertYUnits(yyy[indx+range] - yyy[indx],
1253 
						USER, INCHES, dd);
1254 
			    dxy = hypot(dx, dy);
9279 maechler 1255 
 
13792 pd 1256 
			    /* save the current label for checking overlap */
1257 
			    label2 = allocVector(REALSXP, 8);
9279 maechler 1258 
 
13792 pd 1259 
			    FindCorners(labelDistance, labelHeight, label2,
1260 
					xxx[indx], yyy[indx],
1261 
					xxx[indx+range], yyy[indx+range], dd);
1262 
			    UNPROTECT_PTR(labelList);
1263 
			    labelList = PROTECT(CONS(label2, labelList));
9279 maechler 1264 
 
13792 pd 1265 
			    ddl = FALSE;
1266 
			    /* draw an extra bit of segment if the label
1267 
			       doesn't fill the gap */
7813 murrell 1268 
			    if (closest) {
13792 pd 1269 
				xStart = xxx[indx+range] -
1270 
				    (xxx[indx+range] - xxx[indx]) *
1271 
				    labelDistance / dxy;
1272 
				yStart = yyy[indx+range] -
1273 
				    (yyy[indx+range] - yyy[indx]) *
1274 
				    labelDistance / dxy;
1275 
				if (labelDistance / dxy < 1)
1276 
				    GLine(xxx[indx], yyy[indx],
1277 
					  xStart, yStart,
1278 
					  USER, dd);
7813 murrell 1279 
			    } else {
13792 pd 1280 
				xStart = xxx[indx] +
1281 
				    (xxx[indx+range] - xxx[indx]) *
1282 
				    labelDistance / dxy;
1283 
				yStart = yyy[indx] +
1284 
				    (yyy[indx+range] - yyy[indx]) *
1285 
				    labelDistance / dxy;
1286 
				if (labelDistance / dxy < 1)
1287 
				    GLine(xStart, yStart,
1288 
					  xxx[indx+range], yyy[indx+range],
1289 
					  USER, dd);
7813 murrell 1290 
			    }
13792 pd 1291 
 
1292 
			    /*** Draw contour labels ***/
1293 
			    if (xxx[indx] < xxx[indx+range]) {
1294 
				if (closest) {
1295 
				    ux = xStart;
1296 
				    uy = yStart;
1297 
				    vx = xxx[indx+range];
1298 
				    vy = yyy[indx+range];
1299 
				} else {
1300 
				    ux = xxx[indx];
1301 
				    uy = yyy[indx];
1302 
				    vx = xStart;
1303 
				    vy = yStart;
1304 
				}
7813 murrell 1305 
			    }
13792 pd 1306 
			    else {
1307 
				if (closest) {
1308 
				    ux = xxx[indx+range];
1309 
				    uy = yyy[indx+range];
1310 
				    vx = xStart;
1311 
				    vy = yStart;
1312 
				} else {
1313 
				    ux = xStart;
1314 
				    uy = yStart;
1315 
				    vx = xxx[indx];
1316 
				    vy = yyy[indx];
1317 
				}
1318 
			    }
9279 maechler 1319 
 
13792 pd 1320 
			    if (!ddl) {
7813 murrell 1321 
				/* convert to INCHES for calculation of
9279 maechler 1322 
				   angle to draw text
7813 murrell 1323 
				*/
13792 pd 1324 
				GConvert(&ux, &uy, USER, INCHES, dd);
1325 
				GConvert(&vx, &vy, USER, INCHES, dd);
8062 murrell 1326 
				/* 1, 1 => sans serif, basic font
1327 
				   0, .5 => left, centre justified */
13792 pd 1328 
				if (vectorFonts)
1329 
				    GVText(ux, uy, INCHES, buffer,
1330 
					   typeface, fontindex, 0, .5,
1331 
					   (180 / 3.14) * atan2(vy - uy, vx - ux),
1332 
					   dd);
1333 
				else
1334 
				    GText (ux, uy, INCHES, buffer, 0, .5,
1335 
					   (180 / 3.14) * atan2(vy - uy, vx - ux),
1336 
					   dd);
1337 
			    }
1338 
			} /* if (gotLabel) */
1339 
		    } /* if (method == 0) else ... */
1340 
		} /* if (labelDistance > 0) */
1341 
 
1342 
	    } /* if (drawLabels) */
10886 maechler 1343 
	    else {
1344 
		GPolyline(ns, xxx, yyy, USER, dd);
1345 
	    }
9279 maechler 1346 
 
10886 maechler 1347 
	    GMode(0, dd);
15321 luke 1348 
	    vmaxset(vmax);
10886 maechler 1349 
	} /* while */
1350 
      } /* for(i .. )  for(j ..) */
27811 ripley 1351 
    vmaxset(vmax); /* now we are done with ctr_SegDB */
8872 pd 1352 
    UNPROTECT_PTR(label1); /* pwwwargh! This is messy, but last thing
1353 
			      protected is likely labelList, and that needs
1354 
			      to be preserved across calls */
13792 pd 1355 
}
1284 ihaka 1356 
 
13792 pd 1357 
/* contour(x, y, z, levels, labels, labcex, drawlabels,
1358 
 *         method, vfont, col = col, lty = lty, lwd = lwd)
1359 
 */
1284 ihaka 1360 
SEXP do_contour(SEXP call, SEXP op, SEXP args, SEXP env)
1361 
{
8213 murrell 1362 
    SEXP oargs, c, x, y, z, vfont, col, lty, lwd, labels;
7813 murrell 1363 
    int i, j, nx, ny, nc, ncol, nlty, nlwd;
1364 
    int ltysave, colsave, lwdsave;
1365 
    double cexsave;
1284 ihaka 1366 
    double atom, zmin, zmax;
1367 
    char *vmax, *vmax0;
7813 murrell 1368 
    int method;
10886 maechler 1369 
    Rboolean drawLabels;
1370 
    Rboolean vectorFonts = FALSE;
7813 murrell 1371 
    int typeface = 0;
1372 
    int fontindex = 0;
1373 
    double labcex;
1284 ihaka 1374 
    DevDesc *dd = CurrentDevice();
21061 murrell 1375 
    SEXP result = R_NilValue;
1284 ihaka 1376 
 
1377 
    GCheckState(dd);
1378 
 
1839 ihaka 1379 
    if (length(args) < 4)
5731 ripley 1380 
	errorcall(call, "too few arguments");
1284 ihaka 1381 
 
1382 
    oargs = args;
1383 
 
1384 
    x = CAR(args);
1385 
    internalTypeCheck(call, x, REALSXP);
1386 
    nx = LENGTH(x);
1387 
    args = CDR(args);
1388 
 
1389 
    y = CAR(args);
1390 
    internalTypeCheck(call, y, REALSXP);
1391 
    ny = LENGTH(y);
1392 
    args = CDR(args);
1393 
 
1394 
    z = CAR(args);
1395 
    internalTypeCheck(call, z, REALSXP);
1396 
    args = CDR(args);
1397 
 
13792 pd 1398 
    /* levels */
1284 ihaka 1399 
    c = CAR(args);
1400 
    internalTypeCheck(call, c, REALSXP);
1401 
    nc = LENGTH(c);
1402 
    args = CDR(args);
1403 
 
8213 murrell 1404 
    labels = CAR(args);
1405 
    if (!isNull(labels))
1406 
	internalTypeCheck(call, labels, STRSXP);
1407 
    args = CDR(args);
1408 
 
7813 murrell 1409 
    labcex = asReal(CAR(args));
1410 
    args = CDR(args);
1411 
 
10886 maechler 1412 
    drawLabels = (Rboolean)asLogical(CAR(args));
7813 murrell 1413 
    args = CDR(args);
1414 
 
1415 
    method = asInteger(CAR(args)); args = CDR(args);
1416 
    if (method < 1 || method > 3)
1417 
	errorcall(call, "invalid value for \"method\"");
1418 
 
1419 
    PROTECT(vfont = FixupVFont(CAR(args)));
1420 
    if (!isNull(vfont)) {
10886 maechler 1421 
	vectorFonts = TRUE;
7813 murrell 1422 
	typeface = INTEGER(vfont)[0];
1423 
	fontindex = INTEGER(vfont)[1];
1424 
    }
1425 
    args = CDR(args);
1426 
 
8323 ihaka 1427 
    PROTECT(col = FixupCol(CAR(args), NA_INTEGER));
1284 ihaka 1428 
    ncol = length(col);
8323 ihaka 1429 
    args = CDR(args);
1284 ihaka 1430 
 
17179 murrell 1431 
    PROTECT(lty = FixupLty(CAR(args), Rf_gpptr(dd)->lty));
1284 ihaka 1432 
    nlty = length(lty);
8323 ihaka 1433 
    args = CDR(args);
1284 ihaka 1434 
 
17179 murrell 1435 
    PROTECT(lwd = FixupLwd(CAR(args), Rf_gpptr(dd)->lwd));
7813 murrell 1436 
    nlwd = length(lwd);
8323 ihaka 1437 
    args = CDR(args);
1284 ihaka 1438 
 
1839 ihaka 1439 
    if (nx < 2 || ny < 2)
5731 ripley 1440 
	errorcall(call, "insufficient x or y values");
1284 ihaka 1441 
 
1839 ihaka 1442 
    if (nrows(z) != nx || ncols(z) != ny)
5731 ripley 1443 
	errorcall(call, "dimension mismatch");
1284 ihaka 1444 
 
1839 ihaka 1445 
    if (nc < 1)
5731 ripley 1446 
	errorcall(call, "no contour values");
1284 ihaka 1447 
 
1839 ihaka 1448 
    for (i = 0; i < nx; i++) {
5107 maechler 1449 
	if (!R_FINITE(REAL(x)[i]))
6191 maechler 1450 
	    errorcall(call, "missing x values");
1839 ihaka 1451 
	if (i > 0 && REAL(x)[i] < REAL(x)[i - 1])
5731 ripley 1452 
	    errorcall(call, "increasing x values expected");
1284 ihaka 1453 
    }
1454 
 
1839 ihaka 1455 
    for (i = 0; i < ny; i++) {
5107 maechler 1456 
	if (!R_FINITE(REAL(y)[i]))
6191 maechler 1457 
	    errorcall(call, "missing y values");
1839 ihaka 1458 
	if (i > 0 && REAL(y)[i] < REAL(y)[i - 1])
5731 ripley 1459 
	    errorcall(call, "increasing y values expected");
1284 ihaka 1460 
    }
1461 
 
1839 ihaka 1462 
    for (i = 0; i < nc; i++)
5107 maechler 1463 
	if (!R_FINITE(REAL(c)[i]))
6191 maechler 1464 
	    errorcall(call, "illegal NA contour values");
1284 ihaka 1465 
 
1466 
    zmin = DBL_MAX;
1467 
    zmax = DBL_MIN;
1839 ihaka 1468 
    for (i = 0; i < nx * ny; i++)
5107 maechler 1469 
	if (R_FINITE(REAL(z)[i])) {
1839 ihaka 1470 
	    if (zmax < REAL(z)[i]) zmax =  REAL(z)[i];
1471 
	    if (zmin > REAL(z)[i]) zmin =  REAL(z)[i];
1284 ihaka 1472 
	}
1473 
 
1839 ihaka 1474 
    if (zmin >= zmax) {
1475 
	if (zmin == zmax)
5731 ripley 1476 
	    warning("all z values are equal");
1284 ihaka 1477 
	else
5731 ripley 1478 
	    warning("all z values are NA");
1284 ihaka 1479 
	return R_NilValue;
1480 
    }
1481 
 
13792 pd 1482 
    /* change to 1e-3, reconsidered because of PR#897
1483 
     * but 1e-7, and even  2*DBL_EPSILON do not prevent inf.loop in contour().
1484 
     * maybe something like   16 * DBL_EPSILON * (..).
1485 
     * see also MAX_ns above */
1284 ihaka 1486 
    atom = 1e-3 * (zmax - zmin);
1487 
 
1488 
    /* Initialize the segment data base */
13792 pd 1489 
 
1284 ihaka 1490 
    /* Note we must be careful about resetting */
1491 
    /* the top of the stack, otherwise we run out of */
1492 
    /* memory after a sequence of displaylist replays */
1493 
 
1494 
    vmax0 = vmaxget();
1495 
    ctr_SegDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));
1496 
 
1839 ihaka 1497 
    for (i = 0; i < nx; i++)
1498 
	for (j = 0; j < ny; j++)
1499 
	    ctr_SegDB[i + j * nx] = NULL;
1284 ihaka 1500 
 
1501 
    /* Draw the contours -- note the heap release */
1502 
 
17179 murrell 1503 
    ltysave = Rf_gpptr(dd)->lty;
1504 
    colsave = Rf_gpptr(dd)->col;
1505 
    lwdsave = Rf_gpptr(dd)->lwd;
1506 
    cexsave = Rf_gpptr(dd)->cex;
8872 pd 1507 
    labelList = PROTECT(R_NilValue);
21061 murrell 1508 
 
1509 
 
1510 
    /* draw contour for levels[i] */
3594 ihaka 1511 
    GMode(1, dd);
21061 murrell 1512 
    for (i = 0; i < nc; i++) {
1284 ihaka 1513 
	vmax = vmaxget();
17179 murrell 1514 
	Rf_gpptr(dd)->lty = INTEGER(lty)[i % nlty];
1515 
	if (Rf_gpptr(dd)->lty == NA_INTEGER)
1516 
	    Rf_gpptr(dd)->lty = ltysave;
1517 
	Rf_gpptr(dd)->col = INTEGER(col)[i % ncol];
1518 
	if (Rf_gpptr(dd)->col == NA_INTEGER)
1519 
	    Rf_gpptr(dd)->col = colsave;
1520 
	Rf_gpptr(dd)->lwd = REAL(lwd)[i % nlwd];
1521 
	if (Rf_gpptr(dd)->lwd == NA_REAL)
1522 
	    Rf_gpptr(dd)->lwd = lwdsave;
1523 
	Rf_gpptr(dd)->cex = labcex;
8213 murrell 1524 
	contour(x, nx, y, ny, z, REAL(c)[i], labels, i,
1525 
		drawLabels, method-1,
7813 murrell 1526 
		vectorFonts, typeface, fontindex, atom, dd);
1284 ihaka 1527 
	vmaxset(vmax);
1528 
    }
21061 murrell 1529 
    GMode(0, dd);
1284 ihaka 1530 
    vmaxset(vmax0);
17179 murrell 1531 
    Rf_gpptr(dd)->lty = ltysave;
1532 
    Rf_gpptr(dd)->col = colsave;
1533 
    Rf_gpptr(dd)->lwd = lwdsave;
1534 
    Rf_gpptr(dd)->cex = cexsave;
8872 pd 1535 
    UNPROTECT(5);
1284 ihaka 1536 
    /* NOTE: only record operation if no "error"  */
1537 
    /* NOTE: on replay, call == R_NilValue */
11067 maechler 1538 
    if (GRecording(call))
1284 ihaka 1539 
	recordGraphicOperation(op, oargs, dd);
21061 murrell 1540 
    return result;
1284 ihaka 1541 
}
1542 
 
13792 pd 1543 
 
3866 r 1544 
	/*  F i l l e d   C o n t o u r   P l o t s  */
1545 
 
1546 
	/*  R o s s  I h a k a,  M a r c h  1 9 9 9  */
1547 
 
1548 
static void
1549 
FindCutPoints(double low, double high,
1550 
	      double x1, double y1, double z1,
1551 
	      double x2, double y2, double z2,
1552 
	      double *x, double *y, double *z,
1553 
	      int *npt)
1554 
{
1555 
    double c;
1556 
 
1557 
    if (z1 > z2 ) {
1558 
	if (z2 > high || z1 < low)
1559 
	    return;
1560 
	if (z1 < high) {
1561 
	    x[*npt] = x1;
1562 
	    y[*npt] = y1;
1563 
	    z[*npt] = z1;
1564 
	    ++*npt;
1565 
	}
1566 
	else {
1567 
	    c = (z1 - high) / (z1 - z2);
1568 
	    x[*npt] = x1 + c * (x2 - x1);
1569 
	    y[*npt] = y1;
1570 
	    z[*npt] = z1 + c * (z2 - z1);
1571 
	    ++*npt;
1572 
	}
1573 
	if (z2 > low) {
1574 
	}
1575 
	else {
1576 
	    c = (z2 -low) / (z2 - z1);
1577 
	    x[*npt] = x2 - c * (x2 - x1);
1578 
	    y[*npt] = y1;
1579 
	    z[*npt] = z2 - c * (z2 - z1);
1580 
	    ++*npt;
1581 
	}
1582 
    }
1583 
    else if (z1 < z2) {
1584 
	if (z2 < low || z1 > high)
1585 
	    return;
1586 
	if (z1 > low) {
1587 
	    x[*npt] = x1;
1588 
	    y[*npt] = y1;
1589 
	    z[*npt] = z1;
1590 
	    ++*npt;
1591 
	}
1592 
	else {
1593 
	    c = (z1 - low) / (z1 - z2);
1594 
	    x[*npt] = x1 + c * (x2 - x1);
1595 
	    y[*npt] = y1;
1596 
	    z[*npt] = z1 + c * (z2 - z1);
1597 
	    ++*npt;
1598 
	}
1599 
	if (z2 < high) {
1600 
#ifdef OMIT
1601 
	    /* Don't repeat corner vertices */
1602 
	    x[*npt] = x2;
1603 
	    y[*npt] = y2;
1604 
	    z[*npt] = z2;
1605 
	    ++*npt;
1606 
#endif
1607 
	}
1608 
	else {
1609 
	    c = (z2 - high) / (z2 - z1);
1610 
	    x[*npt] = x2 - c * (x2 - x1);
1611 
	    y[*npt] = y1;
1612 
	    z[*npt] = z2 - c * (z2 - z1);
1613 
	    ++*npt;
1614 
	}
5107 maechler 1615 
    }
3866 r 1616 
    else {
1617 
	if(low <= z1 && z1 <= high) {
1618 
	    x[*npt] = x1;
1619 
	    y[*npt] = y1;
1620 
	    z[*npt] = z1;
1621 
	    ++*npt;
1622 
#ifdef OMIT
1623 
	    /* Don't repeat corner vertices */
1624 
	    x[*npt] = x2;
1625 
	    y[*npt] = y2;
1626 
	    z[*npt] = z2;
1627 
	    ++*npt;
1628 
#endif
5107 maechler 1629 
	}
3866 r 1630 
    }
1631 
}
1632 
 
1633 
/* FIXME - This could pretty easily be adapted to handle NA */
1634 
/* values on the grid.  Just search the diagonals for cutpoints */
1635 
/* instead of the cell sides.  Use the same switch idea as in */
1636 
/* contour above.  There are 5 cases to handle. */
1637 
 
1638 
static void
1639 
FindPolygonVertices(double low, double high,
1640 
		    double x1, double x2, double y1, double y2,
1641 
		    double z11, double z21, double z12, double z22,
1642 
		    double *x, double *y, double *z, int *npt)
1643 
{
1644 
    *npt = 0;
1645 
    FindCutPoints(low, high, x1,  y1,  z11, x2,  y1,  z21, x, y, z, npt);
1646 
    FindCutPoints(low, high, y1,  x2,  z21, y2,  x2,  z22, y, x, z, npt);
1647 
    FindCutPoints(low, high, x2,  y2,  z22, x1,  y2,  z12, x, y, z, npt);
1648 
    FindCutPoints(low, high, y2,  x1,  z12, y1,  x1,  z11, y, x, z, npt);
1649 
}
1650 
 
6098 pd 1651 
/* FIXME: [Code consistency] Use macro for the parallel parts of
1652 
	  do_contour, do_filledcontour & do_image ...
1653 
*/
13792 pd 1654 
 
1655 
/* filledcontour(x, y, z, levels, col) */
3866 r 1656 
SEXP do_filledcontour(SEXP call, SEXP op, SEXP args, SEXP env)
1657 
{
1658 
    SEXP oargs, sx, sy, sz, sc, scol;
1659 
    double *x, *y, *z, *c;
1660 
    unsigned *col;
1661 
    int i, j, k, npt, nx, ny, nz, nc, ncol, colsave, xpdsave;
1662 
    double px[8], py[8], pz[8];
1663 
    DevDesc *dd = CurrentDevice();
1664 
 
1665 
    GCheckState(dd);
1666 
 
1667 
    checkArity(op,args);
1668 
    oargs = args;
1669 
 
1670 
    sx = CAR(args);
1671 
    internalTypeCheck(call, sx, REALSXP);
1672 
    nx = LENGTH(sx);
1673 
    args = CDR(args);
1674 
 
1675 
    sy = CAR(args);
1676 
    internalTypeCheck(call, sy, REALSXP);
1677 
    ny = LENGTH(sy);
1678 
    args = CDR(args);
1679 
 
1680 
    sz = CAR(args);
1681 
    internalTypeCheck(call, sz, REALSXP);
1682 
    nz = length(sz);
1683 
    args = CDR(args);
1684 
 
13792 pd 1685 
    sc = CAR(args);/* levels */
3866 r 1686 
    internalTypeCheck(call, sc, REALSXP);
1687 
    nc = length(sc);
1688 
    args = CDR(args);
1689 
 
19134 ripley 1690 
    if (nx < 2 || ny < 2)
1691 
	errorcall(call, "insufficient x or y values");
1692 
 
1693 
    if (nrows(sz) != nx || ncols(sz) != ny)
1694 
	errorcall(call, "dimension mismatch");
1695 
 
1696 
    if (nc < 1)
1697 
	errorcall(call, "no contour values");
1698 
 
5507 ihaka 1699 
    PROTECT(scol = FixupCol(CAR(args), NA_INTEGER));
3866 r 1700 
    ncol = length(scol);
1701 
 
1702 
    /* Shorthand Pointers */
1703 
 
1704 
    x = REAL(sx);
1705 
    y = REAL(sy);
1706 
    z = REAL(sz);
1707 
    c = REAL(sc);
1708 
    col = (unsigned*)INTEGER(scol);
1709 
 
1710 
    /* Check of grid coordinates */
1711 
    /* We want them to all be finite */
1712 
    /* and in strictly ascending order */
1713 
 
6098 pd 1714 
    if (nx < 1 || ny < 1) goto badxy;
5107 maechler 1715 
    if (!R_FINITE(x[0])) goto badxy;
1716 
    if (!R_FINITE(y[0])) goto badxy;
3866 r 1717 
    for (i = 1; i < nx; i++)
5107 maechler 1718 
	if (!R_FINITE(x[i]) || x[i] <= x[i - 1]) goto badxy;
3866 r 1719 
    for (j = 1; j < ny; j++)
5107 maechler 1720 
	if (!R_FINITE(y[j]) || y[j] <= y[j - 1]) goto badxy;
3866 r 1721 
 
1722 
    /* Check of the contour levels */
1723 
 
5107 maechler 1724 
    if (!R_FINITE(c[0])) goto badlev;
3866 r 1725 
    for (k = 1; k < nc; k++)
5107 maechler 1726 
	if (!R_FINITE(c[k]) || c[k] <= c[k - 1]) goto badlev;
3866 r 1727 
 
17179 murrell 1728 
    colsave = Rf_gpptr(dd)->col;
1729 
    xpdsave = Rf_gpptr(dd)->xpd;
5055 pd 1730 
    /* override par("xpd") and force clipping to plot region */
17179 murrell 1731 
    Rf_gpptr(dd)->xpd = 0;
3866 r 1732 
 
1733 
    GMode(1, dd);
1734 
 
1735 
    for (i = 1; i < nx; i++) {
1736 
	for (j = 1; j < ny; j++) {
1737 
	    for (k = 1; k < nc ; k++) {
1738 
		FindPolygonVertices(c[k - 1], c[k],
1739 
				    x[i - 1], x[i],
1740 
				    y[j - 1], y[j],
1741 
				    z[i - 1 + (j - 1) * nx],
1742 
				    z[i + (j - 1) * nx],
1743 
				    z[i - 1 + j * nx],
1744 
				    z[i + j * nx],
1745 
				    px, py, pz, &npt);
1746 
                if (npt > 2)
1747 
		    GPolygon(npt, px, py, USER, col[(k-1)%ncol],
1748 
			     NA_INTEGER, dd);
1749 
	    }
1750 
	}
1751 
    }
1752 
    GMode(0, dd);
17179 murrell 1753 
    Rf_gpptr(dd)->col = colsave;
1754 
    Rf_gpptr(dd)->xpd = xpdsave;
3866 r 1755 
    R_Visible = 0;
1756 
    UNPROTECT(1);
11067 maechler 1757 
    if (GRecording(call))
3866 r 1758 
	recordGraphicOperation(op, oargs, dd);
1759 
    return R_NilValue;
1760 
 
1761 
 badxy:
6191 maechler 1762 
    errorcall(call, "invalid x / y values or limits");
3866 r 1763 
 badlev:
13792 pd 1764 
    errorcall(call, "invalid contour levels: must be strictly increasing");
5107 maechler 1765 
    return R_NilValue;  /* never used; to keep -Wall happy */
3866 r 1766 
}
1767 
 
1768 
 
1284 ihaka 1769 
	/*  I m a g e   R e n d e r i n g  */
1770 
 
3866 r 1771 
 
8743 ripley 1772 
/* image(x, y, z, col, breaks) */
1284 ihaka 1773 
SEXP do_image(SEXP call, SEXP op, SEXP args, SEXP env)
1774 
{
8743 ripley 1775 
    SEXP oargs, sx, sy, sz, sc;
1776 
    double *x, *y;
1777 
    int *z, tmp;
1284 ihaka 1778 
    unsigned *c;
8743 ripley 1779 
    int i, j, nx, ny, nc, colsave, xpdsave;
1284 ihaka 1780 
    DevDesc *dd = CurrentDevice();
1781 
 
1782 
    GCheckState(dd);
1783 
 
1784 
    checkArity(op,args);
1785 
    oargs = args;
1786 
 
1787 
    sx = CAR(args);
1788 
    internalTypeCheck(call, sx, REALSXP);
1789 
    nx = LENGTH(sx);
1790 
    args = CDR(args);
1791 
 
1792 
    sy = CAR(args);
1793 
    internalTypeCheck(call, sy, REALSXP);
1794 
    ny = LENGTH(sy);
1795 
    args = CDR(args);
1796 
 
1797 
    sz = CAR(args);
8743 ripley 1798 
    internalTypeCheck(call, sz, INTSXP);
1284 ihaka 1799 
    args = CDR(args);
1800 
 
9279 maechler 1801 
    PROTECT(sc = FixupCol(CAR(args), NA_INTEGER));
8743 ripley 1802 
    nc = LENGTH(sc);
1284 ihaka 1803 
 
1804 
    /* Shorthand Pointers */
1805 
 
1806 
    x = REAL(sx);
1807 
    y = REAL(sy);
8743 ripley 1808 
    z = INTEGER(sz);
1284 ihaka 1809 
    c = (unsigned*)INTEGER(sc);
1810 
 
1811 
    /* Check of grid coordinates */
1812 
    /* We want them to all be finite and in strictly ascending order */
1813 
 
6098 pd 1814 
    if (nx < 1 || ny < 1) goto badxy;
5107 maechler 1815 
    if (!R_FINITE(x[0])) goto badxy;
1816 
    if (!R_FINITE(y[0])) goto badxy;
1839 ihaka 1817 
    for (i = 1; i < nx; i++)
5107 maechler 1818 
	if (!R_FINITE(x[i]) || x[i] <= x[i - 1]) goto badxy;
1839 ihaka 1819 
    for (j = 1; j < ny; j++)
5107 maechler 1820 
	if (!R_FINITE(y[j]) || y[j] <= y[j - 1]) goto badxy;
1284 ihaka 1821 
 
17179 murrell 1822 
    colsave = Rf_gpptr(dd)->col;
1823 
    xpdsave = Rf_gpptr(dd)->xpd;
5055 pd 1824 
    /* override par("xpd") and force clipping to plot region */
17179 murrell 1825 
    Rf_gpptr(dd)->xpd = 0;
1284 ihaka 1826 
 
3367 ihaka 1827 
    GMode(1, dd);
9279 maechler 1828 
 
6516 ripley 1829 
    for (i = 0; i < nx - 1 ; i++) {
1830 
	for (j = 0; j < ny - 1; j++) {
8743 ripley 1831 
	    tmp = z[i + j * (nx - 1)];
1832 
	    if (tmp >= 0 && tmp < nc && tmp != NA_INTEGER)
9279 maechler 1833 
		GRect(x[i], y[j], x[i+1], y[j+1], USER, c[tmp],
8743 ripley 1834 
		      NA_INTEGER, dd);
1284 ihaka 1835 
	}
1836 
    }
3367 ihaka 1837 
    GMode(0, dd);
17179 murrell 1838 
    Rf_gpptr(dd)->col = colsave;
1839 
    Rf_gpptr(dd)->xpd = xpdsave;
1284 ihaka 1840 
    R_Visible = 0;
1841 
    UNPROTECT(1);
11067 maechler 1842 
    if (GRecording(call))
1284 ihaka 1843 
	recordGraphicOperation(op, oargs, dd);
1844 
    return R_NilValue;
1845 
 
1846 
  badxy:
6191 maechler 1847 
    errorcall(call, "invalid x / y values or limits");
1284 ihaka 1848 
    return R_NilValue;/* never used; to keep -Wall happy */
1849 
}
1850 
 
1851 
	/*  P e r s p e c t i v e   S u r f a c e   P l o t s  */
1852 
 
1853 
 
3636 r 1854 
/* Set up the light source */
1855 
static double Light[4];
1856 
static double Shade;
19019 ripley 1857 
static Rboolean DoLighting;
3636 r 1858 
 
1859 
static void SetUpLight(double theta, double phi)
1860 
{
3875 hornik 1861 
    double u[4];
3636 r 1862 
    u[0] = 0; u[1] = -1; u[2] = 0; u[3] = 1;
1863 
    SetToIdentity(VT);             /* Initialization */
1864 
    XRotate(-phi);                 /* colatitude rotation */
1865 
    ZRotate(theta);                /* azimuthal rotation */
1866 
    TransVector(u, VT, Light);	   /* transform */
1867 
}
1868 
 
1869 
static double FacetShade(double *u, double *v)
5107 maechler 1870 
{
3636 r 1871 
    double nx, ny, nz, sum;
1872 
    nx = u[1] * v[2] - u[2] * v[1];
1873 
    ny = u[2] * v[0] - u[0] * v[2];
1874 
    nz = u[0] * v[1] - u[1] * v[0];
1875 
    sum = sqrt(nx * nx + ny * ny + nz * nz);
1876 
    if (sum == 0) sum = 1;
1877 
    nx /= sum;
1878 
    ny /= sum;
5107 maechler 1879 
    nz /= sum;
3636 r 1880 
    sum = 0.5 * (nx * Light[0] + ny * Light[1] + nz * Light[2] + 1);
1881 
    return pow(sum, Shade);
1882 
}
1883 
 
1284 ihaka 1884 
 
1592 ihaka 1885 
/* For each facet, determine the farthest point from the eye. */
1284 ihaka 1886 
/* Sorting the facets so that these depths are decreasing */
1592 ihaka 1887 
/* yields an occlusion compatible ordering. */
1888 
/* Note that we ignore z values when doing this. */
1284 ihaka 1889 
 
2342 maechler 1890 
static void DepthOrder(double *z, double *x, double *y, int nx, int ny,
11067 maechler 1891 
		       double *depth, int *indx)
1284 ihaka 1892 
{
1893 
    int i, ii, j, jj, nx1, ny1;
1894 
    Vector3d u, v;
1592 ihaka 1895 
    double d;
1284 ihaka 1896 
    nx1 = nx - 1;
1897 
    ny1 = ny - 1;
1839 ihaka 1898 
    for (i = 0; i < nx1 * ny1; i++)
11067 maechler 1899 
	indx[i] = i;
1839 ihaka 1900 
    for (i = 0; i < nx1; i++)
1901 
	for (j = 0; j < ny1; j++) {
1592 ihaka 1902 
	    d = -DBL_MAX;
1839 ihaka 1903 
	    for (ii = 0; ii <= 1; ii++)
1904 
		for (jj = 0; jj <= 1; jj++) {
1905 
		    u[0] = x[i + ii];
1906 
		    u[1] = y[j + jj];
1592 ihaka 1907 
		    /* Originally I had the following line here: */
1908 
		    /* u[2] = z[i+ii+(j+jj)*nx]; */
1909 
		    /* But this leads to artifacts. */
1910 
		    /* It has been replaced by the following line: */
1911 
		    u[2] = 0;
1284 ihaka 1912 
		    u[3] = 1;
5107 maechler 1913 
		    if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
1592 ihaka 1914 
			TransVector(u, VT, v);
1915 
			if (v[3] > d) d = v[3];
1916 
		    }
1284 ihaka 1917 
		}
9296 ihaka 1918 
	    depth[i+j*nx1] = -d;
2342 maechler 1919 
 
1284 ihaka 1920 
	}
9279 maechler 1921 
    /* Determine the depth ordering of the facets to ensure
1922 
       that they are drawn in an occlusion compatible order. */
11067 maechler 1923 
    rsort_with_index(depth, indx, nx1 * ny1);
1284 ihaka 1924 
}
1925 
 
1592 ihaka 1926 
 
1284 ihaka 1927 
static void DrawFacets(double *z, double *x, double *y, int nx, int ny,
11067 maechler 1928 
		       int *indx, double xs, double ys, double zs,
3636 r 1929 
	               int *col, int ncol, int border)
1284 ihaka 1930 
{
3875 hornik 1931 
    double xx[4], yy[4], shade = 0;
1284 ihaka 1932 
    Vector3d u, v;
1592 ihaka 1933 
    int i, j, k, n, nx1, ny1, icol, nv;
3636 r 1934 
    unsigned int newcol, r, g, b;
1284 ihaka 1935 
    DevDesc *dd;
1936 
    dd = CurrentDevice();
1937 
    nx1 = nx - 1;
1938 
    ny1 = ny - 1;
1939 
    n = nx1 * ny1;
1839 ihaka 1940 
    for (k = 0; k < n; k++) {
1592 ihaka 1941 
	nv = 0;
11067 maechler 1942 
	i = indx[k] % nx1;
1943 
	j = indx[k] / nx1;
1592 ihaka 1944 
	icol = (i + j * nx1) % ncol;
3636 r 1945 
	if (DoLighting) {
1946 
            /* Note we must scale here */
1947 
	    u[0] = xs * (x[i+1] - x[i]);
1948 
	    u[1] = ys * (y[j] - y[j+1]);
1949 
	    u[2] = zs * (z[(i+1)+j*nx] - z[i+(j+1)*nx]);
1950 
	    v[0] = xs * (x[i+1] - x[i]);
1951 
	    v[1] = ys * (y[j+1] - y[j]);
1952 
	    v[2] = zs * (z[(i+1)+(j+1)*nx] - z[i+j*nx]);
1953 
	    shade = FacetShade(u, v);
1954 
	}
1839 ihaka 1955 
	u[0] = x[i]; u[1] = y[j];
1956 
	u[2] = z[i + j * nx]; u[3] = 1;
5107 maechler 1957 
	if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
1592 ihaka 1958 
	    TransVector(u, VT, v);
1959 
	    xx[nv] = v[0] / v[3];
1960 
	    yy[nv] = v[1] / v[3];
1961 
	    nv++;
1962 
	}
1284 ihaka 1963 
 
1839 ihaka 1964 
	u[0] = x[i + 1]; u[1] = y[j];
1965 
	u[2] = z[i + 1 + j * nx]; u[3] = 1;
5107 maechler 1966 
	if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
1592 ihaka 1967 
	    TransVector(u, VT, v);
1968 
	    xx[nv] = v[0] / v[3];
1969 
	    yy[nv] = v[1] / v[3];
1970 
	    nv++;
1971 
	}
1284 ihaka 1972 
 
1839 ihaka 1973 
	u[0] = x[i + 1]; u[1] = y[j + 1];
1974 
	u[2] = z[i + 1 + (j + 1) * nx]; u[3] = 1;
5107 maechler 1975 
	if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
1592 ihaka 1976 
	    TransVector(u, VT, v);
1977 
	    xx[nv] = v[0] / v[3];
1978 
	    yy[nv] = v[1] / v[3];
1979 
	    nv++;
1980 
	}
1284 ihaka 1981 
 
1839 ihaka 1982 
	u[0] = x[i]; u[1] = y[j + 1];
1983 
	u[2] = z[i + (j + 1) * nx]; u[3] = 1;
5107 maechler 1984 
	if (R_FINITE(u[0]) &&  R_FINITE(u[1]) && R_FINITE(u[2])) {
1592 ihaka 1985 
	    TransVector(u, VT, v);
1986 
	    xx[nv] = v[0] / v[3];
1987 
	    yy[nv] = v[1] / v[3];
1988 
	    nv++;
1989 
	}
1284 ihaka 1990 
 
3566 r 1991 
	if (nv > 2) {
3636 r 1992 
	    newcol = col[icol];
1993 
	    if (DoLighting) {
1994 
		r = shade * R_RED(newcol);
1995 
		g = shade * R_GREEN(newcol);
1996 
		b = shade * R_BLUE(newcol);
1997 
		newcol = R_RGB(r, g, b);
1998 
	    }
1999 
	    GPolygon(nv, xx, yy, USER, newcol, border, dd);
3566 r 2000 
	}
1284 ihaka 2001 
    }
2002 
}
2003 
 
1592 ihaka 2004 
 
2342 maechler 2005 
#ifdef NOT_used_currently/*-- out 'def'  (-Wall) --*/
2006 
static void CheckRange(double *x, int n, double min, double max)
1284 ihaka 2007 
{
2008 
    double xmin, xmax;
2009 
    int i;
2010 
    xmin =  DBL_MAX;
2011 
    xmax = -DBL_MAX;
1839 ihaka 2012 
    for (i = 0; i < n; i++)
5107 maechler 2013 
	if (R_FINITE(x[i])) {
1839 ihaka 2014 
	    if (x[i] < xmin) xmin = x[i];
2015 
	    if (x[i] > xmax) xmax = x[i];
1284 ihaka 2016 
	}
1839 ihaka 2017 
    if (xmin < min || xmax > max)
5731 ripley 2018 
	errorcall(gcall, "coordinates outsize specified range");
1284 ihaka 2019 
}
2342 maechler 2020 
#endif
1284 ihaka 2021 
 
2342 maechler 2022 
static void PerspWindow(double *xlim, double *ylim, double *zlim, DevDesc *dd)
1592 ihaka 2023 
{
2024 
    double pin1, pin2, scale, xdelta, ydelta, xscale, yscale, xadd, yadd;
2025 
    double xmax, xmin, ymax, ymin, xx, yy;
2026 
    Vector3d u, v;
2027 
    int i, j, k;
2342 maechler 2028 
 
1592 ihaka 2029 
    xmax = xmin = ymax = ymin = 0;
2030 
    u[3] = 1;
1839 ihaka 2031 
    for (i = 0; i < 2; i++) {
1592 ihaka 2032 
	u[0] = xlim[i];
1839 ihaka 2033 
	for (j = 0; j < 2; j++) {
1592 ihaka 2034 
	    u[1] = ylim[j];
1839 ihaka 2035 
	    for (k = 0; k < 2; k++) {
1592 ihaka 2036 
		u[2] = zlim[k];
2037 
		TransVector(u, VT, v);
2038 
		xx = v[0] / v[3];
2039 
		yy = v[1] / v[3];
1839 ihaka 2040 
		if (xx > xmax) xmax = xx;
2041 
		if (xx < xmin) xmin = xx;
2342 maechler 2042 
		if (yy > ymax) ymax = yy;
2043 
		if (yy < ymin) ymin = yy;
1592 ihaka 2044 
	    }
2045 
	}
2046 
    }
2047 
    pin1 = GConvertXUnits(1.0, NPC, INCHES, dd);
2048 
    pin2 = GConvertYUnits(1.0, NPC, INCHES, dd);
2049 
    xdelta = fabs(xmax - xmin);
2050 
    ydelta = fabs(ymax - ymin);
2051 
    xscale = pin1 / xdelta;
2052 
    yscale = pin2 / ydelta;
2053 
    scale = (xscale < yscale) ? xscale : yscale;
2054 
    xadd = .5 * (pin1 / scale - xdelta);
2055 
    yadd = .5 * (pin2 / scale - ydelta);
2056 
    GScale(xmin - xadd, xmax + xadd, 1, dd);
2057 
    GScale(ymin - yadd, ymax + yadd, 2, dd);
2058 
    GMapWin2Fig(dd);
2059 
}
2060 
 
2061 
static int LimitCheck(double *lim, double *c, double *s)
2062 
{
5107 maechler 2063 
    if (!R_FINITE(lim[0]) || !R_FINITE(lim[1]) || lim[0] >= lim[1])
1592 ihaka 2064 
	return 0;
2065 
    *s = 0.5 * fabs(lim[1] - lim[0]);
2066 
    *c = 0.5 * (lim[1] + lim[0]);
2067 
    return 1;
2068 
}
2069 
 
1671 ihaka 2070 
/* PerspBox: The following code carries out a visibility test */
2071 
/* on the surfaces of the xlim/ylim/zlim box around the plot. */
2072 
/* If front = 0, only the faces with their inside toward the */
2073 
/* eyepoint are drawn.  If front = 1, only the faces with */
2074 
/* their outside toward the eye are drawn.  This lets us carry */
2075 
/* out hidden line removal by drawing any faces which will be */
2076 
/* obscured before the surface, and those which will not be */
2077 
/* obscured after the surface. */
2078 
 
5243 ripley 2079 
/* The vertices of the box */
2080 
static short int Vertex[8][3] = {
3636 r 2081 
    {0, 0, 0},
2082 
    {0, 0, 1},
2083 
    {0, 1, 0},
2084 
    {0, 1, 1},
2085 
    {1, 0, 0},
2086 
    {1, 0, 1},
2087 
    {1, 1, 0},
2088 
    {1, 1, 1},
1647 ihaka 2089 
};
2090 
 
3636 r 2091 
/* The vertices visited when tracing a face */
5243 ripley 2092 
static short int Face[6][4] = {
3636 r 2093 
    {0, 1, 5, 4},
2094 
    {2, 6, 7, 3},
2095 
    {0, 2, 3, 1},
2096 
    {4, 5, 7, 6},
2097 
    {0, 4, 6, 2},
2098 
    {1, 3, 7, 5},
1647 ihaka 2099 
};
2100 
 
3636 r 2101 
/* The edges drawn when tracing a face */
5243 ripley 2102 
static short int Edge[6][4] = {
3636 r 2103 
    { 0, 1, 2, 3},
2104 
    { 4, 5, 6, 7},
2105 
    { 8, 7, 9, 0},
2106 
    { 2,10, 5,11},
2107 
    { 3,11, 4, 8},
2108 
    { 9, 6,10, 1},
2109 
};
2110 
 
2111 
/* Which edges have been drawn previously */
2112 
static char EdgeDone[12];
2113 
 
1647 ihaka 2114 
static void PerspBox(int front, double *x, double *y, double *z, DevDesc *dd)
2115 
{
1671 ihaka 2116 
    Vector3d u0, v0, u1, v1, u2, v2, u3, v3;
2117 
    double d[3], e[3];
6994 pd 2118 
    int f, i, p0, p1, p2, p3, nearby;
17179 murrell 2119 
    int ltysave = Rf_gpptr(dd)->lty;
3636 r 2120 
    if (front)
17179 murrell 2121 
	Rf_gpptr(dd)->lty = LTY_DOTTED;
3636 r 2122 
    else
17179 murrell 2123 
	Rf_gpptr(dd)->lty = LTY_SOLID;
1839 ihaka 2124 
    for (f = 0; f < 6; f++) {
1671 ihaka 2125 
        p0 = Face[f][0];
2126 
        p1 = Face[f][1];
2127 
        p2 = Face[f][2];
2128 
        p3 = Face[f][3];
2129 
 
2130 
	u0[0] = x[Vertex[p0][0]];
2131 
	u0[1] = y[Vertex[p0][1]];
2132 
	u0[2] = z[Vertex[p0][2]];
2133 
	u0[3] = 1;
2134 
	u1[0] = x[Vertex[p1][0]];
2135 
	u1[1] = y[Vertex[p1][1]];
2136 
	u1[2] = z[Vertex[p1][2]];
2137 
	u1[3] = 1;
2138 
	u2[0] = x[Vertex[p2][0]];
2139 
	u2[1] = y[Vertex[p2][1]];
2140 
	u2[2] = z[Vertex[p2][2]];
2141 
	u2[3] = 1;
2142 
	u3[0] = x[Vertex[p3][0]];
2143 
	u3[1] = y[Vertex[p3][1]];
2144 
	u3[2] = z[Vertex[p3][2]];
2145 
	u3[3] = 1;
2342 maechler 2146 
 
1647 ihaka 2147 
	TransVector(u0, VT, v0);
2148 
	TransVector(u1, VT, v1);
1671 ihaka 2149 
	TransVector(u2, VT, v2);
2150 
	TransVector(u3, VT, v3);
2151 
 
3636 r 2152 
	/* Visibility test. */
1839 ihaka 2153 
	/* Determine whether the surface normal is toward the eye. */
3636 r 2154 
	/* Note that we only draw lines once. */
1671 ihaka 2155 
 
1839 ihaka 2156 
        for (i = 0; i < 3; i++) {
1671 ihaka 2157 
	    d[i] = v1[i]/v1[3] - v0[i]/v0[3];
2158 
	    e[i] = v2[i]/v2[3] - v1[i]/v1[3];
2159 
        }
6994 pd 2160 
	nearby = (d[0]*e[1] - d[1]*e[0]) < 0;
1671 ihaka 2161 
 
6994 pd 2162 
	if ((front && nearby) || (!front && !nearby)) {
3636 r 2163 
	    if (!EdgeDone[Edge[f][0]]++)
2164 
		GLine(v0[0]/v0[3], v0[1]/v0[3],
2165 
		      v1[0]/v1[3], v1[1]/v1[3], USER, dd);
2166 
	    if (!EdgeDone[Edge[f][1]]++)
2167 
		GLine(v1[0]/v1[3], v1[1]/v1[3],
2168 
		      v2[0]/v2[3], v2[1]/v2[3], USER, dd);
2169 
	    if (!EdgeDone[Edge[f][2]]++)
2170 
		GLine(v2[0]/v2[3], v2[1]/v2[3],
2171 
		      v3[0]/v3[3], v3[1]/v3[3], USER, dd);
2172 
	    if (!EdgeDone[Edge[f][3]]++)
2173 
		GLine(v3[0]/v3[3], v3[1]/v3[3],
2174 
		      v0[0]/v0[3], v0[1]/v0[3], USER, dd);
1647 ihaka 2175 
	}
2176 
    }
17179 murrell 2177 
    Rf_gpptr(dd)->lty = ltysave;
1647 ihaka 2178 
}
2179 
 
7813 murrell 2180 
/* PerspAxes:
2181 
 */
2182 
 
2183 
/* Starting vertex for possible axes */
2184 
static short int AxisStart[8] = { 0, 0, 2, 4, 0, 4, 2, 6 };
2185 
 
2186 
/* Tick vector for possible axes */
2187 
static short int TickVector[8][3] = {
2188 
    {0, -1, -1},
2189 
    {-1, 0, -1},
2190 
    {0, 1, -1},
2191 
    {1, 0, -1},
2192 
    {-1, -1, 0},
2193 
    {1, -1, 0},
2194 
    {-1, 1, 0},
2195 
    {1, 1, 0}};
2196 
 
8374 ripley 2197 
static int lowest(double y1, double y2, double y3, double y4) {
7813 murrell 2198 
    return ((y1 <= y2) && (y1 <= y3) && (y1 <= y4));
2199 
}
2200 
 
8374 ripley 2201 
static double labelAngle(double x1, double y1, double x2, double y2) {
8036 murrell 2202 
    double dx, dy;
2203 
    double angle;
2204 
    dx = fabs(x2 - x1);
9279 maechler 2205 
    if (x2 > x1)
8036 murrell 2206 
	dy = y2 - y1;
2207 
    else
2208 
	dy = y1 - y2;
2209 
    if (dx == 0) {
2210 
	if (dy > 0)
2211 
	    angle = 90;
2212 
	else
2213 
	    angle = 270;
2214 
    } else {
2215 
	angle = (180 / M_PI) * atan2(dy, dx);
2216 
    }
2217 
    return angle;
2218 
}
2219 
 
9279 maechler 2220 
static void PerspAxis(double *x, double *y, double *z,
2221 
		      int axis, int axisType, int nTicks, int tickType,
8036 murrell 2222 
		      char *label, DevDesc *dd) {
8331 hornik 2223 
    Vector3d u1, u2, u3, v1, v2, v3;
8062 murrell 2224 
    double tickLength = .03; /* proportion of axis length */
12256 pd 2225 
    double min, max, d_frac;
8578 murrell 2226 
    double *range = NULL; /* -Wall */
7813 murrell 2227 
    double axp[3];
7824 ripley 2228 
    int nint, i;
7813 murrell 2229 
    SEXP at, lab;
2230 
    switch (axisType) {
2231 
    case 0:
12256 pd 2232 
	min = x[0];	max = x[1];	range = x;	break;
7813 murrell 2233 
    case 1:
12256 pd 2234 
	min = y[0];	max = y[1];	range = y;	break;
7813 murrell 2235 
    case 2:
12256 pd 2236 
	min = z[0];	max = z[1];	range = z;	break;
7813 murrell 2237 
    }
12256 pd 2238 
    d_frac = 0.1*(max - min);
2239 
    nint = nTicks - 1; if(!nint) nint++;
2240 
    i = nint;
7813 murrell 2241 
    GPretty(&min, &max, &nint);
12256 pd 2242 
    /* GPretty() rarely gives values too much outside range ..
2243 
       2D axis() clip these, we play cheaper */
2244 
    while((min < range[0] - d_frac || range[1] + d_frac < max) && i < 20) {
2245 
	nint = ++i;
2246 
	min = range[0];
2247 
	max = range[1];
2248 
	GPretty(&min, &max, &nint);
2249 
    }
7813 murrell 2250 
    axp[0] = min;
2251 
    axp[1] = max;
2252 
    axp[2] = nint;
8036 murrell 2253 
    /* Do the following calculations for both ticktypes */
2254 
    switch (axisType) {
2255 
    case 0:
9279 maechler 2256 
	u1[0] = min;
8036 murrell 2257 
	u1[1] = y[Vertex[AxisStart[axis]][1]];
2258 
	u1[2] = z[Vertex[AxisStart[axis]][2]];
2259 
	break;
2260 
    case 1:
2261 
	u1[0] = x[Vertex[AxisStart[axis]][0]];
2262 
	u1[1] = min;
2263 
	u1[2] = z[Vertex[AxisStart[axis]][2]];
2264 
	break;
2265 
    case 2:
2266 
	u1[0] = x[Vertex[AxisStart[axis]][0]];
2267 
	u1[1] = y[Vertex[AxisStart[axis]][1]];
2268 
	u1[2] = min;
2269 
	break;
2270 
    }
2271 
    u1[0] = u1[0] + tickLength*(x[1]-x[0])*TickVector[axis][0];
2272 
    u1[1] = u1[1] + tickLength*(y[1]-y[0])*TickVector[axis][1];
2273 
    u1[2] = u1[2] + tickLength*(z[1]-z[0])*TickVector[axis][2];
2274 
    u1[3] = 1;
2275 
    switch (axisType) {
2276 
    case 0:
9279 maechler 2277 
	u2[0] = max;
8036 murrell 2278 
	u2[1] = u1[1];
2279 
	u2[2] = u1[2];
2280 
	break;
2281 
    case 1:
2282 
	u2[0] = u1[0];
2283 
	u2[1] = max;
2284 
	u2[2] = u1[2];
2285 
	break;
2286 
    case 2:
2287 
	u2[0] = u1[0];
2288 
	u2[1] = u1[1];
2289 
	u2[2] = max;
2290 
	break;
2291 
    }
2292 
    u2[3] = 1;
2293 
    /* The axis label has to be further out for "detailed" ticks
2294 
       in order to leave room for the tick labels */
2295 
    switch (tickType) {
9279 maechler 2296 
    case 1: /* "simple": just an arrow parallel to axis, indicating direction
8036 murrell 2297 
	       of increase */
2298 
	u3[0] = u1[0] + tickLength*(x[1]-x[0])*TickVector[axis][0];
2299 
	u3[1] = u1[1] + tickLength*(y[1]-y[0])*TickVector[axis][1];
2300 
	u3[2] = u1[2] + tickLength*(z[1]-z[0])*TickVector[axis][2];
2301 
	break;
2302 
    case 2:
2303 
	u3[0] = u1[0] + 2.5*tickLength*(x[1]-x[0])*TickVector[axis][0];
2304 
	u3[1] = u1[1] + 2.5*tickLength*(y[1]-y[0])*TickVector[axis][1];
2305 
	u3[2] = u1[2] + 2.5*tickLength*(z[1]-z[0])*TickVector[axis][2];
2306 
	break;
2307 
    }
2308 
    switch (axisType) {
2309 
    case 0:
9279 maechler 2310 
	u3[0] = (min + max)/2;
8036 murrell 2311 
	break;
2312 
    case 1:
2313 
	u3[1] = (min + max)/2;
2314 
	break;
2315 
    case 2:
2316 
	u3[2] = (min + max)/2;
2317 
	break;
2318 
    }
2319 
    u3[3] = 1;
2320 
    TransVector(u1, VT, v1);
2321 
    TransVector(u2, VT, v2);
2322 
    TransVector(u3, VT, v3);
2323 
    /* Draw axis label */
9279 maechler 2324 
    GText(v3[0]/v3[3], v3[1]/v3[3], USER, label, .5, .5,
8036 murrell 2325 
	  labelAngle(v1[0]/v1[3], v1[1]/v1[3], v2[0]/v2[3], v2[1]/v2[3]),
2326 
	  dd);
2327 
    /* Draw axis ticks */
2328 
    switch (tickType) {
9279 maechler 2329 
    case 1: /* "simple": just an arrow parallel to axis, indicating direction
8036 murrell 2330 
	       of increase */
2331 
	/* arrow head is 0.25 inches long, with angle 30 degrees,
2332 
	   and drawn at v2 end of line */
2333 
	GArrow(v1[0]/v1[3], v1[1]/v1[3],
9279 maechler 2334 
	       v2[0]/v2[3], v2[1]/v2[3], USER,
8036 murrell 2335 
	       0.1, 10, 2, dd);
2336 
	break;
2337 
    case 2: /* "detailed": normal ticks as per 2D plots */
10886 maechler 2338 
	PROTECT(at = CreateAtVector(axp, range, 7, FALSE));
8036 murrell 2339 
	PROTECT(lab = labelformat(at));
9279 maechler 2340 
	for (i=0; i<length(at); i++) {
8036 murrell 2341 
	    switch (axisType) {
2342 
	    case 0:
9279 maechler 2343 
		u1[0] = REAL(at)[i];
8036 murrell 2344 
		u1[1] = y[Vertex[AxisStart[axis]][1]];
2345 
		u1[2] = z[Vertex[AxisStart[axis]][2]];
2346 
		break;
2347 
	    case 1:
2348 
		u1[0] = x[Vertex[AxisStart[axis]][0]];
2349 
		u1[1] = REAL(at)[i];
2350 
		u1[2] = z[Vertex[AxisStart[axis]][2]];
2351 
		break;
2352 
	    case 2:
2353 
		u1[0] = x[Vertex[AxisStart[axis]][0]];
2354 
		u1[1] = y[Vertex[AxisStart[axis]][1]];
2355 
		u1[2] = REAL(at)[i];
2356 
		break;
2357 
	    }
2358 
	    u1[3] = 1;
2359 
	    u2[0] = u1[0] + tickLength*(x[1]-x[0])*TickVector[axis][0];
2360 
	    u2[1] = u1[1] + tickLength*(y[1]-y[0])*TickVector[axis][1];
2361 
	    u2[2] = u1[2] + tickLength*(z[1]-z[0])*TickVector[axis][2];
2362 
	    u2[3] = 1;
2363 
	    u3[0] = u2[0] + tickLength*(x[1]-x[0])*TickVector[axis][0];
2364 
	    u3[1] = u2[1] + tickLength*(y[1]-y[0])*TickVector[axis][1];
2365 
	    u3[2] = u2[2] + tickLength*(z[1]-z[0])*TickVector[axis][2];
2366 
	    u3[3] = 1;
2367 
	    TransVector(u1, VT, v1);
2368 
	    TransVector(u2, VT, v2);
2369 
	    TransVector(u3, VT, v3);
2370 
	    /* Draw tick line */
2371 
	    GLine(v1[0]/v1[3], v1[1]/v1[3],
2372 
		  v2[0]/v2[3], v2[1]/v2[3], USER, dd);
2373 
	    /* Draw tick label */
10172 luke 2374 
	    GText(v3[0]/v3[3], v3[1]/v3[3], USER, CHAR(STRING_ELT(lab, i)),
8036 murrell 2375 
		  .5, .5, 0, dd);
7813 murrell 2376 
	}
8036 murrell 2377 
	UNPROTECT(2);
2378 
	break;
7813 murrell 2379 
    }
2380 
}
2381 
 
18912 murrell 2382 
/* Determine the transformed (x, y) coordinates (in USER space)
2383 
 * for the four corners of the x-y plane of the persp plot
2384 
 * These will be used to determine which sides of the persp
2385 
 * plot to label with axes
2386 
 * The strategy is to determine which corner has the lowest y-value
2387 
 * to decide which of the x- and y-axes to label AND which corner
2388 
 * has the lowest x-value to decide which of the z-axes to label
2389 
 */
9279 maechler 2390 
static void PerspAxes(double *x, double *y, double *z,
8036 murrell 2391 
                      char *xlab, char *ylab, char *zlab,
2392 
		      int nTicks, int tickType, DevDesc *dd) {
7824 ripley 2393 
    int xAxis=0, yAxis=0, zAxis=0; /* -Wall */
7813 murrell 2394 
    int xpdsave;
18912 murrell 2395 
    Vector3d u0, u1, u2, u3;
7813 murrell 2396 
    Vector3d v0, v1, v2, v3;
18912 murrell 2397 
    u0[0] = x[0];
2398 
    u0[1] = y[0];
2399 
    u0[2] = z[0];
2400 
    u0[3] = 1;
2401 
    u1[0] = x[1];
2402 
    u1[1] = y[0];
2403 
    u1[2] = z[0];
2404 
    u1[3] = 1;
2405 
    u2[0] = x[0];
2406 
    u2[1] = y[1];
2407 
    u2[2] = z[0];
2408 
    u2[3] = 1;
2409 
    u3[0] = x[1];
2410 
    u3[1] = y[1];
2411 
    u3[2] = z[0];
2412 
    u3[3] = 1;
7813 murrell 2413 
    TransVector(u0, VT, v0);
2414 
    TransVector(u1, VT, v1);
2415 
    TransVector(u2, VT, v2);
2416 
    TransVector(u3, VT, v3);
2417 
 
2418 
    /* to fit in the axis labels */
17179 murrell 2419 
    xpdsave = Rf_gpptr(dd)->xpd;
2420 
    Rf_gpptr(dd)->xpd = 1;
7813 murrell 2421 
 
2422 
    /* Figure out which X and Y axis to draw */
2423 
    if (lowest(v0[1]/v0[3], v1[1]/v1[3], v2[1]/v2[3], v3[1]/v3[3])) {
2424 
	xAxis = 0;
2425 
	yAxis = 1;
2426 
    } else if (lowest(v1[1]/v1[3], v0[1]/v0[3], v2[1]/v2[3], v3[1]/v3[3])) {
2427 
	xAxis = 0;
2428 
	yAxis = 3;
2429 
    } else if (lowest(v2[1]/v2[3], v1[1]/v1[3], v0[1]/v0[3], v3[1]/v3[3])) {
2430 
	xAxis = 2;
2431 
	yAxis = 1;
2432 
    } else if (lowest(v3[1]/v3[3], v1[1]/v1[3], v2[1]/v2[3], v0[1]/v0[3])) {
2433 
	xAxis = 2;
2434 
	yAxis = 3;
9279 maechler 2435 
    } else
7813 murrell 2436 
	warning("Axis orientation not calculated");
8036 murrell 2437 
    PerspAxis(x, y, z, xAxis, 0, nTicks, tickType, xlab, dd);
2438 
    PerspAxis(x, y, z, yAxis, 1, nTicks, tickType, ylab, dd);
7813 murrell 2439 
    /* Figure out which Z axis to draw */
2440 
    if (lowest(v0[0]/v0[3], v1[0]/v1[3], v2[0]/v2[3], v3[0]/v3[3])) {
2441 
	zAxis = 4;
2442 
    } else if (lowest(v1[0]/v1[3], v0[0]/v0[3], v2[0]/v2[3], v3[0]/v3[3])) {
2443 
	zAxis = 5;
2444 
    } else if (lowest(v2[0]/v2[3], v1[0]/v1[3], v0[0]/v0[3], v3[0]/v3[3])) {
2445 
	zAxis = 6;
2446 
    } else if (lowest(v3[0]/v3[3], v1[0]/v1[3], v2[0]/v2[3], v0[0]/v0[3])) {
2447 
	zAxis = 7;
2448 
    } else
2449 
	warning("Axes orientation not calculated");
8036 murrell 2450 
    PerspAxis(x, y, z, zAxis, 2, nTicks, tickType, zlab, dd);
9279 maechler 2451 
 
17179 murrell 2452 
    Rf_gpptr(dd)->xpd = xpdsave;
9279 maechler 2453 
}
2454 
 
1284 ihaka 2455 
SEXP do_persp(SEXP call, SEXP op, SEXP args, SEXP env)
2456 
{
1858 ihaka 2457 
    SEXP x, y, z, xlim, ylim, zlim;
11067 maechler 2458 
    SEXP depth, indx, originalArgs;
8036 murrell 2459 
    SEXP col, border, xlab, ylab, zlab;
3636 r 2460 
    double theta, phi, r, d;
2461 
    double ltheta, lphi;
2462 
    double expand, xc, yc, zc, xs, ys, zs;
8036 murrell 2463 
    int i, j, scale, ncol, dobox, doaxes, nTicks, tickType;
5513 ihaka 2464 
    DevDesc *dd;
1284 ihaka 2465 
 
12256 pd 2466 
    if (length(args) < 24)
5731 ripley 2467 
	errorcall(call, "too few parameters");
1284 ihaka 2468 
    gcall = call;
2469 
    originalArgs = args;
2470 
 
1592 ihaka 2471 
    PROTECT(x = coerceVector(CAR(args), REALSXP));
5731 ripley 2472 
    if (length(x) < 2) errorcall(call, "invalid x argument");
1592 ihaka 2473 
    args = CDR(args);
2474 
 
2475 
    PROTECT(y = coerceVector(CAR(args), REALSXP));
5731 ripley 2476 
    if (length(y) < 2) errorcall(call, "invalid y argument");
1592 ihaka 2477 
    args = CDR(args);
2478 
 
1284 ihaka 2479 
    PROTECT(z = coerceVector(CAR(args), REALSXP));
1592 ihaka 2480 
    if (!isMatrix(z) || nrows(z) != length(x) || ncols(z) != length(y))
5731 ripley 2481 
	errorcall(call, "invalid z argument");
1284 ihaka 2482 
    args = CDR(args);
2483 
 
1592 ihaka 2484 
    PROTECT(xlim = coerceVector(CAR(args), REALSXP));
5731 ripley 2485 
    if (length(xlim) != 2) errorcall(call, "invalid xlim argument");
1284 ihaka 2486 
    args = CDR(args);
2487 
 
1592 ihaka 2488 
    PROTECT(ylim = coerceVector(CAR(args), REALSXP));
5731 ripley 2489 
    if (length(ylim) != 2) errorcall(call, "invalid ylim argument");
1284 ihaka 2490 
    args = CDR(args);
2491 
 
1592 ihaka 2492 
    PROTECT(zlim = coerceVector(CAR(args), REALSXP));
5731 ripley 2493 
    if (length(zlim) != 2) errorcall(call, "invalid zlim argument");
1592 ihaka 2494 
    args = CDR(args);
1284 ihaka 2495 
 
1592 ihaka 2496 
    /* Checks on x/y/z Limits */
1284 ihaka 2497 
 
1592 ihaka 2498 
    if (!LimitCheck(REAL(xlim), &xc, &xs))
5731 ripley 2499 
	errorcall(call, "invalid x limits");
1592 ihaka 2500 
    if (!LimitCheck(REAL(ylim), &yc, &ys))
5731 ripley 2501 
	errorcall(call, "invalid y limits");
1592 ihaka 2502 
    if (!LimitCheck(REAL(zlim), &zc, &zs))
5731 ripley 2503 
	errorcall(call, "invalid z limits");
1592 ihaka 2504 
 
6415 maechler 2505 
    theta = asReal(CAR(args));	args = CDR(args);
2506 
    phi	  = asReal(CAR(args));	args = CDR(args);
2507 
    r	= asReal(CAR(args));	args = CDR(args);
2508 
    d	= asReal(CAR(args));	args = CDR(args);
2509 
    scale  = asLogical(CAR(args)); args = CDR(args);
2510 
    expand = asReal(CAR(args)); args = CDR(args);
2511 
    col	   = CAR(args);		args = CDR(args);
2512 
    border = CAR(args);		args = CDR(args);
2513 
    ltheta = asReal(CAR(args)); args = CDR(args);
2514 
    lphi   = asReal(CAR(args)); args = CDR(args);
2515 
    Shade  = asReal(CAR(args)); args = CDR(args);
2516 
    dobox  = asLogical(CAR(args)); args = CDR(args);
7813 murrell 2517 
    doaxes = asLogical(CAR(args)); args = CDR(args);
8036 murrell 2518 
    nTicks = asInteger(CAR(args)); args = CDR(args);
2519 
    tickType = asInteger(CAR(args)); args = CDR(args);
2520 
    xlab = CAR(args); args = CDR(args);
2521 
    ylab = CAR(args); args = CDR(args);
2522 
    zlab = CAR(args); args = CDR(args);
19135 ripley 2523 
    if (!isString(xlab) || length(xlab) < 1)
2524 
	error("`xlab' must be a character vector of length 1");
2525 
    if (!isString(ylab) || length(ylab) < 1)
2526 
	error("`ylab' must be a character vector of length 1");
2527 
    if (!isString(zlab) || length(zlab) < 1)
2528 
	error("`zlab' must be a character vector of length 1");
1592 ihaka 2529 
 
5107 maechler 2530 
    if (R_FINITE(Shade) && Shade <= 0) Shade = 1;
2531 
    if (R_FINITE(ltheta) && R_FINITE(lphi) && R_FINITE(Shade))
19019 ripley 2532 
	DoLighting = TRUE;
3636 r 2533 
    else
19019 ripley 2534 
	DoLighting = FALSE;
3636 r 2535 
 
1592 ihaka 2536 
    if (!scale) {
2537 
	double s;
2538 
	s = xs;
2539 
	if (s < ys) s = ys;
2540 
	if (s < zs) s = zs;
2541 
	xs = s; ys = s; zs = s;
2542 
    }
2543 
 
2544 
    /* Parameter Checks */
2545 
 
5107 maechler 2546 
    if (!R_FINITE(theta) || !R_FINITE(phi) || !R_FINITE(r) || !R_FINITE(d) ||
1592 ihaka 2547 
	d < 0 || r < 0)
6191 maechler 2548 
	errorcall(call, "invalid viewing parameters");
5107 maechler 2549 
    if (!R_FINITE(expand) || expand < 0)
6191 maechler 2550 
	errorcall(call, "invalid expand value");
1592 ihaka 2551 
    if (scale == NA_LOGICAL)
2552 
	scale = 0;
9279 maechler 2553 
    if ((nTicks == NA_INTEGER) || (nTicks < 0))
8036 murrell 2554 
	errorcall(call, "invalid nticks value");
9279 maechler 2555 
    if ((tickType == NA_INTEGER) || (tickType < 1) || (tickType > 2))
8036 murrell 2556 
	errorcall(call, "invalid ticktype value");
1592 ihaka 2557 
 
11112 maechler 2558 
    dd = GNewPlot(GRecording(call));
5513 ihaka 2559 
 
17179 murrell 2560 
    PROTECT(col = FixupCol(col, Rf_gpptr(dd)->bg));
5513 ihaka 2561 
    ncol = LENGTH(col);
6191 maechler 2562 
    if (ncol < 1) errorcall(call, "invalid col specification");
19019 ripley 2563 
    if(!R_OPAQUE(INTEGER(col)[0])) DoLighting = FALSE;
17179 murrell 2564 
    PROTECT(border = FixupCol(border, Rf_gpptr(dd)->fg));
6191 maechler 2565 
    if (length(border) < 1) errorcall(call, "invalid border specification");
5513 ihaka 2566 
 
1284 ihaka 2567 
    GSetState(1, dd);
2568 
    GSavePars(dd);
19912 duncan 2569 
    ProcessInlinePars(args, dd, call);
1839 ihaka 2570 
    if (length(border) > 1)
17179 murrell 2571 
	Rf_gpptr(dd)->fg = INTEGER(border)[0];
2572 
    Rf_gpptr(dd)->xlog = Rf_gpptr(dd)->ylog = FALSE;
1284 ihaka 2573 
 
3636 r 2574 
    /* Set up the light vector (if any) */
2575 
    if (DoLighting)
2576 
        SetUpLight(ltheta, lphi);
2577 
 
2578 
    /* Mark box edges as undrawn */
2579 
    for (i = 0; i< 12; i++)
2580 
	EdgeDone[i] = 0;
2581 
 
1284 ihaka 2582 
    /* Specify the viewing transformation. */
2342 maechler 2583 
 
1592 ihaka 2584 
    SetToIdentity(VT);             /* Initialization */
2585 
    Translate(-xc, -yc, -zc);      /* center at the origin */
2586 
    Scale(1/xs, 1/ys, expand/zs);  /* scale extents to [-1,1] */
2587 
    XRotate(-90.0);                /* rotate x-y plane to horizontal */
2588 
    YRotate(-theta);               /* azimuthal rotation */
2589 
    XRotate(phi);                  /* elevation rotation */
2590 
    Translate(0.0, 0.0, -r - d);   /* translate the eyepoint to the origin */
2591 
    Perspective(d);                /* perspective */
2342 maechler 2592 
 
1592 ihaka 2593 
    /* Specify the plotting window. */
2594 
    /* Here we map the vertices of the cube */
2595 
    /* [xmin,xmax]*[ymin,ymax]*[zmin,zmax] */
2596 
    /* to the screen and then chose a window */
2597 
    /* which is symmetric about (0,0). */
2598 
 
2599 
    PerspWindow(REAL(xlim), REAL(ylim), REAL(zlim), dd);
2600 
 
2601 
    /* Compute facet order. */
2602 
 
1284 ihaka 2603 
    PROTECT(depth = allocVector(REALSXP, (nrows(z) - 1)*(ncols(z) - 1)));
11067 maechler 2604 
    PROTECT(indx = allocVector(INTSXP, (nrows(z) - 1)*(ncols(z) - 1)));
1284 ihaka 2605 
    DepthOrder(REAL(z), REAL(x), REAL(y), nrows(z), ncols(z),
11067 maechler 2606 
	       REAL(depth), INTEGER(indx));
1284 ihaka 2607 
 
12256 pd 2608 
    /* Now we order the facets by depth and then draw them back to front.
2609 
     * This is the "painters" algorithm. */
1284 ihaka 2610 
 
3594 ihaka 2611 
    GMode(1, dd);
9279 maechler 2612 
 
7813 murrell 2613 
    if (dobox) {
3636 r 2614 
        PerspBox(0, REAL(xlim), REAL(ylim), REAL(zlim), dd);
20377 ripley 2615 
	if (doaxes) {
2616 
	    SEXP xl = STRING_ELT(xlab, 0), yl = STRING_ELT(ylab, 0),
2617 
		zl = STRING_ELT(zlab, 0);
9279 maechler 2618 
	    PerspAxes(REAL(xlim), REAL(ylim), REAL(zlim),
20377 ripley 2619 
		      (xl == NA_STRING)? "" : CHAR(xl),
2620 
		      (yl == NA_STRING)? "" : CHAR(yl),
2621 
		      (zl == NA_STRING)? "" : CHAR(zl),
8036 murrell 2622 
		      nTicks, tickType, dd);
20377 ripley 2623 
	}
7813 murrell 2624 
    }
1647 ihaka 2625 
 
11067 maechler 2626 
    DrawFacets(REAL(z), REAL(x), REAL(y), nrows(z), ncols(z), INTEGER(indx),
3636 r 2627 
	       1/xs, 1/ys, expand/zs,
2628 
	       INTEGER(col), ncol, INTEGER(border)[0]);
1284 ihaka 2629 
 
3636 r 2630 
    if (dobox)
2631 
        PerspBox(1, REAL(xlim), REAL(ylim), REAL(zlim), dd);
3594 ihaka 2632 
    GMode(0, dd);
1647 ihaka 2633 
 
1284 ihaka 2634 
    GRestorePars(dd);
3636 r 2635 
    UNPROTECT(10);
11067 maechler 2636 
    if (GRecording(call))
1284 ihaka 2637 
        recordGraphicOperation(op, originalArgs, dd);
1647 ihaka 2638 
 
2639 
    PROTECT(x = allocVector(REALSXP, 16));
2640 
    PROTECT(y = allocVector(INTSXP, 2));
1839 ihaka 2641 
    for (i = 0; i < 4; i++)
12256 pd 2642 
	for (j = 0; j < 4; j++)
2643 
	    REAL(x)[i + j * 4] = VT[i][j];
1647 ihaka 2644 
    INTEGER(y)[0] = 4;
2645 
    INTEGER(y)[1] = 4;
2646 
    setAttrib(x, R_DimSymbol, y);
2647 
    UNPROTECT(2);
1284 ihaka 2648 
    return x;
2649 
}
7813 murrell 2650