The R Project SVN R

Rev

Rev 59258 | Rev 60473 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 59258 Rev 59267
Line 33... Line 33...
33 
 
 33 
 
34 
/* filled contours and perspective plots were originally here,
 34 
/* filled contours and perspective plots were originally here,
35 
   now in graphics/src/plot3d.c .
 35 
   now in graphics/src/plot3d.c .
36 
 */
 36 
 */
37 
 
 37 
 
38 
/* Stuff for labels on contour plots
 - 
 
39 
   Originally written by Nicholas Hildreth
 - 
 
40 
   Adapted by Paul Murrell
 - 
 
41 
*/
 - 
 
42 
 
 - 
 
43 
static
 - 
 
44 
void FindCorners(double width, double height, SEXP label,
 - 
 
45 
		 double x0, double y0, double x1, double y1,
 - 
 
46 
		 pGEDevDesc dd) {
 - 
 
47 
    double delta = height / width;
 - 
 
48 
    double dx = GConvertXUnits(x1 - x0, USER, INCHES, dd) * delta;
 - 
 
49 
    double dy = GConvertYUnits(y1 - y0, USER, INCHES, dd) * delta;
 - 
 
50 
    dx = GConvertYUnits(dx, INCHES, USER, dd);
 - 
 
51 
    dy = GConvertXUnits(dy, INCHES, USER, dd);
 - 
 
52 
 
 - 
 
53 
    REAL(label)[0] = x0 + dy;
 - 
 
54 
    REAL(label)[4] = y0 - dx;
 - 
 
55 
    REAL(label)[1] = x0 - dy;
 - 
 
56 
    REAL(label)[5] = y0 + dx;
 - 
 
57 
    REAL(label)[3] = x1 + dy;
 - 
 
58 
    REAL(label)[7] = y1 - dx;
 - 
 
59 
    REAL(label)[2] = x1 - dy;
 - 
 
60 
    REAL(label)[6] = y1 + dx;
 - 
 
61 
}
 - 
 
62 
 
 - 
 
63 
static
 - 
 
64 
int TestLabelIntersection(SEXP label1, SEXP label2) {
 - 
 
65 
 
 - 
 
66 
    int i, j, l1, l2;
 - 
 
67 
    double Ax, Bx, Ay, By, ax, ay, bx, by;
 - 
 
68 
    double dom;
 - 
 
69 
    double result1, result2;
 - 
 
70 
 
 - 
 
71 
    for (i = 0; i < 4; i++) {
 - 
 
72 
	Ax = REAL(label1)[i];
 - 
 
73 
	Ay = REAL(label1)[i+4];
 - 
 
74 
	Bx = REAL(label1)[(i+1)%4];
 - 
 
75 
	By = REAL(label1)[(i+1)%4+4];
 - 
 
76 
	for (j = 0; j < 4; j++) {
 - 
 
77 
	    ax = REAL(label2)[j];
 - 
 
78 
	    ay = REAL(label2)[j+4];
 - 
 
79 
	    bx = REAL(label2)[(j+1)%4];
 - 
 
80 
	    by = REAL(label2)[(j+1)%4+4];
 - 
 
81 
 
 - 
 
82 
	    dom = Bx*by - Bx*ay - Ax*by + Ax*ay - bx*By + bx*Ay + ax*By - ax*Ay;
 - 
 
83 
	    if (dom == 0.0) {
 - 
 
84 
		result1 = -1;
 - 
 
85 
		result2 = -1;
 - 
 
86 
	    }
 - 
 
87 
	    else {
 - 
 
88 
		result1 = (bx*Ay - ax*Ay - ay*bx - Ax*by + Ax*ay + by*ax) / dom;
 - 
 
89 
 
 - 
 
90 
		if (bx - ax == 0.0) {
 - 
 
91 
		    if (by - ay == 0.0)
 - 
 
92 
			result2 = -1;
 - 
 
93 
		    else
 - 
 
94 
			result2 = (Ay + (By - Ay) * result1 - ay) / (by - ay);
 - 
 
95 
		}
 - 
 
96 
		else
 - 
 
97 
		    result2 = (Ax + (Bx - Ax) * result1 - ax) / (bx - ax);
 - 
 
98 
 
 - 
 
99 
	    }
 - 
 
100 
	    l1 = (result1 >= 0.0) && (result1 <= 1.0);
 - 
 
101 
	    l2 = (result2 >= 0.0) && (result2 <= 1.0);
 - 
 
102 
	    if (l1 && l2) return 1;
 - 
 
103 
	}
 - 
 
104 
    }
 - 
 
105 
 
 - 
 
106 
    return 0;
 - 
 
107 
}
 - 
 
108 
 
 - 
 
109 
/*** Checks whether a label window is inside view region ***/
 - 
 
110 
static int LabelInsideWindow(SEXP label, pGEDevDesc dd) {
 - 
 
111 
    int i = 0;
 - 
 
112 
    double x, y;
 - 
 
113 
 
 - 
 
114 
    while (i < 4) {
 - 
 
115 
	x = REAL(label)[i];
 - 
 
116 
	y = REAL(label)[i+4];
 - 
 
117 
	GConvert(&x, &y, USER, NDC, dd);
 - 
 
118 
	/*	x = GConvertXUnits(REAL(label)[i], USER, NDC, dd);
 - 
 
119 
		y = GConvertYUnits(REAL(label)[i+4], USER, NDC, dd); */
 - 
 
120 
 
 - 
 
121 
	if ((x < 0) || (x > 1) ||
 - 
 
122 
	    (y < 0) || (y > 1))
 - 
 
123 
	    return 1;
 - 
 
124 
	i += 1;
 - 
 
125 
    }
 - 
 
126 
    return 0;
 - 
 
127 
}
 - 
 
128 
 
 - 
 
129 
 
 - 
 
130 
	/*  C o n t o u r   P l o t t i n g  */
 - 
 
131 
 
 - 
 
132 
typedef struct SEG {
 - 
 
133 
    struct SEG *next;
 - 
 
134 
    double x0;
 - 
 
135 
    double y0;
 - 
 
136 
    double x1;
 - 
 
137 
    double y1;
 - 
 
138 
} SEG, *SEGP;
 - 
 
139 
 
 - 
 
140 
static SEGP *ctr_SegDB;
 - 
 
141 
 
 - 
 
142 
static int ctr_intersect(double z0, double z1, double zc, double *f)
 - 
 
143 
{
 - 
 
144 
    if ((z0 - zc) * (z1 - zc) < 0.0) {
 - 
 
145 
	*f = (zc - z0) / (z1 -	z0);
 - 
 
146 
	return 1;
 - 
 
147 
    }
 - 
 
148 
    return 0;
 - 
 
149 
}
 - 
 
150 
 
 - 
 
151 
static SEGP ctr_newseg(double x0, double y0, double x1, double y1, SEGP prev)
 - 
 
152 
{
 - 
 
153 
    SEGP seg = (SEGP)R_alloc(1, sizeof(SEG));
 - 
 
154 
    seg->x0 = x0;
 - 
 
155 
    seg->y0 = y0;
 - 
 
156 
    seg->x1 = x1;
 - 
 
157 
    seg->y1 = y1;
 - 
 
158 
    seg->next = prev;
 - 
 
159 
    return seg;
 - 
 
160 
}
 - 
 
161 
 
 - 
 
162 
static void ctr_swapseg(SEGP seg)
 - 
 
163 
{
 - 
 
164 
    double x, y;
 - 
 
165 
    x = seg->x0;
 - 
 
166 
    y = seg->y0;
 - 
 
167 
    seg->x0 = seg->x1;
 - 
 
168 
    seg->y0 = seg->y1;
 - 
 
169 
    seg->x1 = x;
 - 
 
170 
    seg->y1 = y;
 - 
 
171 
}
 - 
 
172 
 
 - 
 
173 
	/* ctr_segdir(): Determine the entry direction to the next cell */
 - 
 
174 
	/* and update the cell indices */
 - 
 
175 
 
 - 
 
176 
#define XMATCH(x0,x1) (fabs(x0-x1) == 0)
 - 
 
177 
#define YMATCH(y0,y1) (fabs(y0-y1) == 0)
 - 
 
178 
 
 - 
 
179 
static int ctr_segdir(double xend, double yend, double *x, double *y,
 - 
 
180 
		      int *i, int *j, int nx, int ny)
 - 
 
181 
{
 - 
 
182 
    if (YMATCH(yend, y[*j])) {
 - 
 
183 
	if (*j == 0)
 - 
 
184 
	    return 0;
 - 
 
185 
	*j = *j - 1;
 - 
 
186 
	return 3;
 - 
 
187 
    }
 - 
 
188 
    if (XMATCH(xend, x[*i])) {
 - 
 
189 
	if (*i == 0)
 - 
 
190 
	    return 0;
 - 
 
191 
	*i = *i - 1;
 - 
 
192 
	return 4;
 - 
 
193 
    }
 - 
 
194 
    if (YMATCH(yend, y[*j + 1])) {
 - 
 
195 
	if (*j >= ny - 1)
 - 
 
196 
	    return 0;
 - 
 
197 
	*j = *j + 1;
 - 
 
198 
	return 1;
 - 
 
199 
    }
 - 
 
200 
    if (XMATCH(xend, x[*i + 1])) {
 - 
 
201 
	if (*i >= nx - 1)
 - 
 
202 
	    return 0;
 - 
 
203 
	*i = *i + 1;
 - 
 
204 
	return 2;
 - 
 
205 
    }
 - 
 
206 
    return 0;
 - 
 
207 
}
 - 
 
208 
 
 - 
 
209 
/* Search seglist for a segment with endpoint (xend, yend). */
 - 
 
210 
/* The cell entry direction is dir, and if tail=1/0 we are */
 - 
 
211 
/* building the tail/head of a contour.	 The matching segment */
 - 
 
212 
/* is pointed to by seg and the updated segment list (with */
 - 
 
213 
/* the matched segment stripped) is returned by the funtion. */
 - 
 
214 
 
 - 
 
215 
static SEGP ctr_segupdate(double xend, double yend, int dir, Rboolean tail,
 - 
 
216 
			  SEGP seglist, SEGP* seg)
 - 
 
217 
{
 - 
 
218 
    if (seglist == NULL) {
 - 
 
219 
	*seg = NULL;
 - 
 
220 
	return NULL;
 - 
 
221 
    }
 - 
 
222 
    switch (dir) {
 - 
 
223 
    case 1:
 - 
 
224 
    case 3:
 - 
 
225 
	if (YMATCH(yend,seglist->y0)) {
 - 
 
226 
	    if (!tail)
 - 
 
227 
		ctr_swapseg(seglist);
 - 
 
228 
	    *seg = seglist;
 - 
 
229 
	    return seglist->next;
 - 
 
230 
	}
 - 
 
231 
	if (YMATCH(yend,seglist->y1)) {
 - 
 
232 
	    if (tail)
 - 
 
233 
		ctr_swapseg(seglist);
 - 
 
234 
	    *seg = seglist;
 - 
 
235 
	    return seglist->next;
 - 
 
236 
	}
 - 
 
237 
	break;
 - 
 
238 
    case 2:
 - 
 
239 
    case 4:
 - 
 
240 
	if (XMATCH(xend,seglist->x0)) {
 - 
 
241 
	    if (!tail)
 - 
 
242 
		ctr_swapseg(seglist);
 - 
 
243 
	    *seg = seglist;
 - 
 
244 
	    return seglist->next;
 - 
 
245 
	}
 - 
 
246 
	if (XMATCH(xend,seglist->x1)) {
 - 
 
247 
	    if (tail)
 - 
 
248 
		ctr_swapseg(seglist);
 - 
 
249 
	    *seg = seglist;
 - 
 
250 
	    return seglist->next;
 - 
 
251 
	}
 - 
 
252 
	break;
 - 
 
253 
    }
 - 
 
254 
    seglist->next = ctr_segupdate(xend, yend, dir, tail, seglist->next, seg);
 - 
 
255 
    return seglist;
 - 
 
256 
}
 - 
 
257 
 
 - 
 
258 
/* labelList, label1, and label2 are all SEXPs rather than being allocated
 - 
 
259 
   using R_alloc because they need to persist across calls to contour().
 - 
 
260 
   In do_contour() there is a vmaxget() ... vmaxset() around each call to
 - 
 
261 
   contour() to release all of the memory used in the drawing of the
 - 
 
262 
   contour _lines_ at each contour level.  We need to keep track of the
 - 
 
263 
   contour _labels_ for _all_ contour levels, hence we have to use a
 - 
 
264 
   different memory allocation mechanism.
 - 
 
265 
*/
 - 
 
266 
static SEXP labelList;
 - 
 
267 
 
 - 
 
268 
static
 - 
 
269 
double distFromEdge(double *xxx, double *yyy, int iii, pGEDevDesc dd) {
 - 
 
270 
    return fmin2(fmin2(xxx[iii]-gpptr(dd)->usr[0], gpptr(dd)->usr[1]-xxx[iii]),
 - 
 
271 
		 fmin2(yyy[iii]-gpptr(dd)->usr[2], gpptr(dd)->usr[3]-yyy[iii]));
 - 
 
272 
}
 - 
 
273 
 
 - 
 
274 
static
 - 
 
275 
Rboolean useStart(double *xxx, double *yyy, int ns, pGEDevDesc dd) {
 - 
 
276 
    if (distFromEdge(xxx, yyy, 0, dd) < distFromEdge(xxx, yyy, ns-1, dd))
 - 
 
277 
	return TRUE;
 - 
 
278 
    else
 - 
 
279 
	return FALSE;
 - 
 
280 
}
 - 
 
281 
 
 - 
 
282 
static
 - 
 
283 
int findGapUp(double *xxx, double *yyy, int ns, double labelDistance,
 - 
 
284 
	      pGEDevDesc dd) {
 - 
 
285 
    double dX, dY;
 - 
 
286 
    double dXC, dYC;
 - 
 
287 
    double distanceSum = 0;
 - 
 
288 
    int n = 0;
 - 
 
289 
    int jjj = 1;
 - 
 
290 
    while ((jjj < ns) && (distanceSum < labelDistance)) {
 - 
 
291 
	/* Find a gap big enough for the label
 - 
 
292 
	   use several segments if necessary
 - 
 
293 
	*/
 - 
 
294 
	dX = xxx[jjj] - xxx[jjj - n - 1]; /* jjj - n - 1 == 0 */
 - 
 
295 
	dY = yyy[jjj] - yyy[jjj - n - 1];
 - 
 
296 
	dXC = GConvertXUnits(dX, USER, INCHES, dd);
 - 
 
297 
	dYC = GConvertYUnits(dY, USER, INCHES, dd);
 - 
 
298 
	distanceSum = hypot(dXC, dYC);
 - 
 
299 
	jjj++;
 - 
 
300 
	n++;
 - 
 
301 
    }
 - 
 
302 
    if (distanceSum < labelDistance)
 - 
 
303 
	return 0;
 - 
 
304 
    else
 - 
 
305 
	return n;
 - 
 
306 
}
 - 
 
307 
 
 - 
 
308 
static
 - 
 
309 
int findGapDown(double *xxx, double *yyy, int ns, double labelDistance,
 - 
 
310 
		pGEDevDesc dd) {
 - 
 
311 
    double dX, dY;
 - 
 
312 
    double dXC, dYC;
 - 
 
313 
    double distanceSum = 0;
 - 
 
314 
    int n = 0;
 - 
 
315 
    int jjj = ns - 2;
 - 
 
316 
    while ((jjj > -1) && (distanceSum < labelDistance)) {
 - 
 
317 
	/* Find a gap big enough for the label
 - 
 
318 
	   use several segments if necessary
 - 
 
319 
	*/
 - 
 
320 
	dX = xxx[jjj] - xxx[jjj + n + 1]; /*jjj + n + 1 == ns -1 */
 - 
 
321 
	dY = yyy[jjj] - yyy[jjj + n + 1];
 - 
 
322 
	dXC = GConvertXUnits(dX, USER, INCHES, dd);
 - 
 
323 
	dYC = GConvertYUnits(dY, USER, INCHES, dd);
 - 
 
324 
	distanceSum = hypot(dXC, dYC);
 - 
 
325 
	jjj--;
 - 
 
326 
	n++;
 - 
 
327 
    }
 - 
 
328 
    if (distanceSum < labelDistance)
 - 
 
329 
	return 0;
 - 
 
330 
    else
 - 
 
331 
	return n;
 - 
 
332 
}
 - 
 
333 
 
 - 
 
334 
/*
 - 
 
335 
 * Generate a list of segments for a single level
 - 
 
336 
 */
 - 
 
337 
static SEGP* contourLines(double *x, int nx, double *y, int ny,
 - 
 
338 
			 double *z, double zc, double atom)
 - 
 
339 
{
 - 
 
340 
    double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];
 - 
 
341 
    int i, j, k, l, m, nacode;
 - 
 
342 
    SEGP seglist;
 - 
 
343 
    SEGP *segmentDB;
 - 
 
344 
    /* Initialize the segment data base */
 - 
 
345 
    /* Note we must be careful about resetting */
 - 
 
346 
    /* the top of the stack, otherwise we run out of */
 - 
 
347 
    /* memory after a sequence of displaylist replays */
 - 
 
348 
    /*
 - 
 
349 
     * This reset is done out in GEcontourLines
 - 
 
350 
     */
 - 
 
351 
    segmentDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));
 - 
 
352 
    for (i = 0; i < nx; i++)
 - 
 
353 
	for (j = 0; j < ny; j++)
 - 
 
354 
	    segmentDB[i + j * nx] = NULL;
 - 
 
355 
    for (i = 0; i < nx - 1; i++) {
 - 
 
356 
	xl = x[i];
 - 
 
357 
	xh = x[i + 1];
 - 
 
358 
	for (j = 0; j < ny - 1; j++) {
 - 
 
359 
	    yl = y[j];
 - 
 
360 
	    yh = y[j + 1];
 - 
 
361 
	    k = i + j * nx;
 - 
 
362 
	    zll = z[k];
 - 
 
363 
	    zhl = z[k + 1];
 - 
 
364 
	    zlh = z[k + nx];
 - 
 
365 
	    zhh = z[k + nx + 1];
 - 
 
366 
 
 - 
 
367 
	    /* If the value at a corner is exactly equal to a contour level,
 - 
 
368 
	     * change that value by a tiny amount */
 - 
 
369 
 
 - 
 
370 
	    if (zll == zc) zll += atom;
 - 
 
371 
	    if (zhl == zc) zhl += atom;
 - 
 
372 
	    if (zlh == zc) zlh += atom;
 - 
 
373 
	    if (zhh == zc) zhh += atom;
 - 
 
374 
#ifdef DEBUG_contour
 38 
#include "contour-common.h"
375 
	    /* Haven't seen this happening (MM): */
 - 
 
376 
	    if (zll == zc) REprintf(" [%d,%d] ll: %g\n",i,j, zll);
 - 
 
377 
	    if (zhl == zc) REprintf(" [%d,%d] hl: %g\n",i,j, zhl);
 - 
 
378 
	    if (zlh == zc) REprintf(" [%d,%d] lh: %g\n",i,j, zlh);
 - 
 
379 
	    if (zhh == zc) REprintf(" [%d,%d] hh: %g\n",i,j, zhh);
 - 
 
380 
#endif
 - 
 
381 
	    /* Check for intersections with sides */
 - 
 
382 
 
 - 
 
383 
	    nacode = 0;
 - 
 
384 
	    if (R_FINITE(zll)) nacode += 1;
 - 
 
385 
	    if (R_FINITE(zhl)) nacode += 2;
 - 
 
386 
	    if (R_FINITE(zlh)) nacode += 4;
 - 
 
387 
	    if (R_FINITE(zhh)) nacode += 8;
 - 
 
388 
 
 - 
 
389 
	    k = 0;
 - 
 
390 
	    switch (nacode) {
 - 
 
391 
	    case 15:
 - 
 
392 
		if (ctr_intersect(zll, zhl, zc, &f)) {
 - 
 
393 
		    xx[k] = xl + f * (xh - xl);
 - 
 
394 
		    yy[k] = yl; k++;
 - 
 
395 
		}
 - 
 
396 
		if (ctr_intersect(zll, zlh, zc, &f)) {
 - 
 
397 
		    yy[k] = yl + f * (yh - yl);
 - 
 
398 
		    xx[k] = xl; k++;
 - 
 
399 
		}
 - 
 
400 
		if (ctr_intersect(zhl, zhh, zc, &f)) {
 - 
 
401 
		    yy[k] = yl + f * (yh - yl);
 - 
 
402 
		    xx[k] = xh; k++;
 - 
 
403 
		}
 - 
 
404 
		if (ctr_intersect(zlh, zhh, zc, &f)) {
 - 
 
405 
		    xx[k] = xl + f * (xh - xl);
 - 
 
406 
		    yy[k] = yh; k++;
 - 
 
407 
		}
 - 
 
408 
		break;
 - 
 
409 
	    case 14:
 - 
 
410 
		if (ctr_intersect(zhl, zhh, zc, &f)) {
 - 
 
411 
		    yy[k] = yl + f * (yh - yl);
 - 
 
412 
		    xx[k] = xh; k++;
 - 
 
413 
		}
 - 
 
414 
		if (ctr_intersect(zlh, zhh, zc, &f)) {
 - 
 
415 
		    xx[k] = xl + f * (xh - xl);
 - 
 
416 
		    yy[k] = yh; k++;
 - 
 
417 
		}
 - 
 
418 
		if (ctr_intersect(zlh, zhl, zc, &f)) {
 - 
 
419 
		    xx[k] = xl + f * (xh - xl);
 - 
 
420 
		    yy[k] = yh + f * (yl - yh);
 - 
 
421 
		    k++;
 - 
 
422 
		}
 - 
 
423 
		break;
 - 
 
424 
	    case 13:
 - 
 
425 
		if (ctr_intersect(zll, zlh, zc, &f)) {
 - 
 
426 
		    yy[k] = yl + f * (yh - yl);
 - 
 
427 
		    xx[k] = xl; k++;
 - 
 
428 
		}
 - 
 
429 
		if (ctr_intersect(zlh, zhh, zc, &f)) {
 - 
 
430 
		    xx[k] = xl + f * (xh - xl);
 - 
 
431 
		    yy[k] = yh; k++;
 - 
 
432 
		}
 - 
 
433 
		if (ctr_intersect(zll, zhh, zc, &f)) {
 - 
 
434 
		    xx[k] = xl + f * (xh - xl);
 - 
 
435 
		    yy[k] = yl + f * (yh - yl);
 - 
 
436 
		    k++;
 - 
 
437 
		}
 - 
 
438 
		break;
 - 
 
439 
	    case 11:
 - 
 
440 
		if (ctr_intersect(zhl, zhh, zc, &f)) {
 - 
 
441 
		    yy[k] = yl + f * (yh - yl);
 - 
 
442 
		    xx[k] = xh; k++;
 - 
 
443 
		}
 - 
 
444 
		if (ctr_intersect(zll, zhl, zc, &f)) {
 - 
 
445 
		    xx[k] = xl + f * (xh - xl);
 - 
 
446 
		    yy[k] = yl; k++;
 - 
 
447 
		}
 - 
 
448 
		if (ctr_intersect(zll, zhh, zc, &f)) {
 - 
 
449 
		    xx[k] = xl + f * (xh - xl);
 - 
 
450 
		    yy[k] = yl + f * (yh - yl);
 - 
 
451 
		    k++;
 - 
 
452 
		}
 - 
 
453 
		break;
 - 
 
454 
	    case 7:
 - 
 
455 
		if (ctr_intersect(zll, zlh, zc, &f)) {
 - 
 
456 
		    yy[k] = yl + f * (yh - yl);
 - 
 
457 
		    xx[k] = xl; k++;
 - 
 
458 
		}
 - 
 
459 
		if (ctr_intersect(zll, zhl, zc, &f)) {
 - 
 
460 
		    xx[k] = xl + f * (xh - xl);
 - 
 
461 
		    yy[k] = yl; k++;
 - 
 
462 
		}
 - 
 
463 
		if (ctr_intersect(zlh, zhl, zc, &f)) {
 - 
 
464 
		    xx[k] = xl + f * (xh - xl);
 - 
 
465 
		    yy[k] = yh + f * (yl - yh);
 - 
 
466 
		    k++;
 - 
 
467 
		}
 - 
 
468 
		break;
 - 
 
469 
	    }
 - 
 
470 
 
 - 
 
471 
	    /* We now have k(=2,4) endpoints */
 - 
 
472 
	    /* Decide which to join */
 - 
 
473 
 
 - 
 
474 
	    seglist = NULL;
 - 
 
475 
 
 - 
 
476 
	    if (k > 0) {
 - 
 
477 
		if (k == 2) {
 - 
 
478 
		    seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
 - 
 
479 
		}
 - 
 
480 
		else if (k == 4) {
 - 
 
481 
		    for (k = 3; k >= 1; k--) {
 - 
 
482 
			m = k;
 - 
 
483 
			xl = xx[k];
 - 
 
484 
			for (l = 0; l < k; l++) {
 - 
 
485 
			    if (xx[l] > xl) {
 - 
 
486 
				xl = xx[l];
 - 
 
487 
				m = l;
 - 
 
488 
			    }
 - 
 
489 
			}
 - 
 
490 
			if (m != k) {
 - 
 
491 
			    xl = xx[k];
 - 
 
492 
			    yl = yy[k];
 - 
 
493 
			    xx[k] = xx[m];
 - 
 
494 
			    yy[k] = yy[m];
 - 
 
495 
			    xx[m] = xl;
 - 
 
496 
			    yy[m] = yl;
 - 
 
497 
			}
 - 
 
498 
		    }
 - 
 
499 
		    seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
 - 
 
500 
		    seglist = ctr_newseg(xx[2], yy[2], xx[3], yy[3], seglist);
 - 
 
501 
		}
 - 
 
502 
		else error("k != 2 or 4");
 - 
 
503 
	    }
 - 
 
504 
	    segmentDB[i + j * nx] = seglist;
 - 
 
505 
	}
 - 
 
506 
    }
 - 
 
507 
    return segmentDB;
 - 
 
508 
}
 - 
 
509 
 
 39 
 
510 
#define CONTOUR_LIST_STEP 100
 40 
#define CONTOUR_LIST_STEP 100
511 
#define CONTOUR_LIST_LEVEL 0
 41 
#define CONTOUR_LIST_LEVEL 0
512 
#define CONTOUR_LIST_X 1
 42 
#define CONTOUR_LIST_X 1
513 
#define CONTOUR_LIST_Y 2
 43 
#define CONTOUR_LIST_Y 2
Line 728... Line 258...
728 
	mainlist = VECTOR_ELT(container, 0);
 258 
	mainlist = VECTOR_ELT(container, 0);
729 
    UNPROTECT(1);  /* UNPROTECT container */
 259 
    UNPROTECT(1);  /* UNPROTECT container */
730 
    return mainlist;
 260 
    return mainlist;
731 
}
 261 
}
732 
 
 262 
 
- 
 
 263 
/* This is for contourLines() in package grDevices */
733 
SEXP attribute_hidden do_contourLines(SEXP call, SEXP op, SEXP args, SEXP env)
 264 
SEXP attribute_hidden do_contourLines(SEXP call, SEXP op, SEXP args, SEXP env)
734 
{
 265 
{
735 
    SEXP c, x, y, z;
 266 
    SEXP c, x, y, z;
736 
    int nx, ny, nc;
 267 
    int nx, ny, nc;
737 
 
 268 
 
Line 755... Line 286...
755 
    nc = LENGTH(c);
 286 
    nc = LENGTH(c);
756 
    args = CDR(args);
 287 
    args = CDR(args);
757 
 
 288 
 
758 
    return GEcontourLines(REAL(x), nx, REAL(y), ny, REAL(z), REAL(c), nc);
 289 
    return GEcontourLines(REAL(x), nx, REAL(y), ny, REAL(z), REAL(c), nc);
759 
}
 290 
}
760 
 
 - 
 
761 
/*
 - 
 
762 
 * The *base graphics* function contour() and the *general base*
 - 
 
763 
 * function contourLines() use the same code to generate contour lines
 - 
 
764 
 * (i.e., the function contourLines())
 - 
 
765 
 *
 - 
 
766 
 * I had a look at extracting the code that draws the labels
 - 
 
767 
 * into a *general base* function
 - 
 
768 
 * (e.g., into some sort of labelLines() function),
 - 
 
769 
 * but the code is too base-graphics-specific (e.g., one of the
 - 
 
770 
 * labelling methods seeks the location closest to the edge of the
 - 
 
771 
 * plotting region) so I've left it alone for now.
 - 
 
772 
 *
 - 
 
773 
 * This does mean that the contourLines() function is part of the
 - 
 
774 
 * graphics engine, but the contour() function is part of the
 - 
 
775 
 * base graphics system.
 - 
 
776 
 */
 - 
 
777 
 
 - 
 
778 
static void contour(SEXP x, int nx, SEXP y, int ny, SEXP z,
 - 
 
779 
		    double zc,
 - 
 
780 
		    SEXP labels, int cnum,
 - 
 
781 
		    Rboolean drawLabels, int method,
 - 
 
782 
		    double atom, pGEDevDesc dd)
 - 
 
783 
{
 - 
 
784 
/* draw a contour for one given contour level 'zc' */
 - 
 
785 
 
 - 
 
786 
    const void *vmax;
 - 
 
787 
 
 - 
 
788 
    double xend, yend;
 - 
 
789 
    int i, ii, j, jj, ns, dir;
 - 
 
790 
    SEGP seglist, seg, s, start, end;
 - 
 
791 
    double *xxx, *yyy;
 - 
 
792 
 
 - 
 
793 
    double variance, dX, dY, deltaX, deltaY;
 - 
 
794 
    double dXC, dYC;
 - 
 
795 
    int range=0, indx=0, n; /* -Wall */
 - 
 
796 
    double lowestVariance;
 - 
 
797 
    double squareSum;
 - 
 
798 
    int iii, jjj;
 - 
 
799 
    double distanceSum, labelDistance, avgGradient;
 - 
 
800 
    char buffer[255];
 - 
 
801 
    int result;
 - 
 
802 
    double ux, uy, vx, vy;
 - 
 
803 
    double xStart, yStart;
 - 
 
804 
    double dx, dy, dxy;
 - 
 
805 
    double labelHeight;
 - 
 
806 
    SEXP label1 = PROTECT(allocVector(REALSXP, 8));
 - 
 
807 
    SEXP label2;
 - 
 
808 
    SEXP lab;
 - 
 
809 
    Rboolean gotLabel = FALSE;
 - 
 
810 
    Rboolean ddl;/* Don't draw label -- currently unused, i.e. always FALSE*/
 - 
 
811 
 
 - 
 
812 
#ifdef DEBUG_contour
 - 
 
813 
    Rprintf("contour(lev = %g):\n", zc);
 - 
 
814 
#endif
 - 
 
815 
 
 - 
 
816 
    vmax = vmaxget();
 - 
 
817 
    ctr_SegDB = contourLines(REAL(x), nx, REAL(y), ny, REAL(z), zc, atom);
 - 
 
818 
    /* we need to keep ctr_SegDB available, so vmaxset(vmax); was wrong */
 - 
 
819 
 
 - 
 
820 
    /* The segment database is now assembled. */
 - 
 
821 
    /* Begin following contours. */
 - 
 
822 
    /* 1. Grab a segment */
 - 
 
823 
    /* 2. Follow its tail */
 - 
 
824 
    /* 3. Follow its head */
 - 
 
825 
    /* 4. Draw the contour */
 - 
 
826 
 
 - 
 
827 
    for (i = 0; i < nx - 1; i++)
 - 
 
828 
      for (j = 0; j < ny - 1; j++) {
 - 
 
829 
	while ((seglist = ctr_SegDB[i + j * nx])) {
 - 
 
830 
	    ii = i; jj = j;
 - 
 
831 
	    start = end = seglist;
 - 
 
832 
	    ctr_SegDB[i + j * nx] = seglist->next;
 - 
 
833 
	    xend = seglist->x1;
 - 
 
834 
	    yend = seglist->y1;
 - 
 
835 
	    while ((dir = ctr_segdir(xend, yend, REAL(x), REAL(y),
 - 
 
836 
				     &ii, &jj, nx, ny))) {
 - 
 
837 
		ctr_SegDB[ii + jj * nx]
 - 
 
838 
		    = ctr_segupdate(xend, yend, dir, TRUE,/* = tail */
 - 
 
839 
				    ctr_SegDB[ii + jj * nx], &seg);
 - 
 
840 
		if (!seg) break;
 - 
 
841 
		end->next = seg;
 - 
 
842 
		end = seg;
 - 
 
843 
		xend = end->x1;
 - 
 
844 
		yend = end->y1;
 - 
 
845 
	    }
 - 
 
846 
	    end->next = NULL; /* <<< new for 1.2.3 */
 - 
 
847 
	    ii = i; jj = j;
 - 
 
848 
	    xend = seglist->x0;
 - 
 
849 
	    yend = seglist->y0;
 - 
 
850 
	    while ((dir = ctr_segdir(xend, yend, REAL(x), REAL(y),
 - 
 
851 
				     &ii, &jj, nx, ny))) {
 - 
 
852 
		ctr_SegDB[ii + jj * nx]
 - 
 
853 
		    = ctr_segupdate(xend, yend, dir, FALSE,/* ie. head */
 - 
 
854 
				    ctr_SegDB[ii+jj*nx], &seg);
 - 
 
855 
		if (!seg) break;
 - 
 
856 
		seg->next = start;
 - 
 
857 
		start = seg;
 - 
 
858 
		xend = start->x0;
 - 
 
859 
		yend = start->y0;
 - 
 
860 
	    }
 - 
 
861 
 
 - 
 
862 
	    /* ns := #{segments of polyline} -- need to allocate */
 - 
 
863 
	    s = start;
 - 
 
864 
	    ns = 0;
 - 
 
865 
	    /* max_contour_segments: prevent inf.loop (shouldn't be needed) */
 - 
 
866 
	    while (s && ns < max_contour_segments) {
 - 
 
867 
		ns++;
 - 
 
868 
		s = s->next;
 - 
 
869 
	    }
 - 
 
870 
	    if(ns == max_contour_segments)
 - 
 
871 
		warning(_("contour(): circular/long seglist -- set %s > %d?"),
- 
 
872 
		        "options(\"max.contour.segments\")", max_contour_segments);
 - 
 
873 
 
 - 
 
874 
	    /* contour midpoint : use for labelling sometime (not yet!)
 - 
 
875 
	       int ns2;
 - 
 
876 
	       if (ns > 3) ns2 = ns/2; else ns2 = -1;
 - 
 
877 
	    */
 - 
 
878 
 
 - 
 
879 
	    vmax = vmaxget();
 - 
 
880 
	    xxx = (double *) R_alloc(ns + 1, sizeof(double));
 - 
 
881 
	    yyy = (double *) R_alloc(ns + 1, sizeof(double));
 - 
 
882 
	    /* now have the space, go through again: */
 - 
 
883 
	    s = start;
 - 
 
884 
	    ns = 0;
 - 
 
885 
	    xxx[ns] = s->x0;
 - 
 
886 
	    yyy[ns++] = s->y0;
 - 
 
887 
	    while (s->next && ns < max_contour_segments) {
 - 
 
888 
		s = s->next;
 - 
 
889 
		xxx[ns] = s->x0;
 - 
 
890 
		yyy[ns++] = s->y0;
 - 
 
891 
	    }
 - 
 
892 
	    xxx[ns] = s->x1;
 - 
 
893 
	    yyy[ns++] = s->y1;
 - 
 
894 
#ifdef DEBUG_contour
 - 
 
895 
	    Rprintf("  [%2d,%2d]: (x,y)[1:%d] = ", i,j, ns);
 - 
 
896 
	    if(ns >= 5)
 - 
 
897 
		Rprintf(" (%g,%g), (%g,%g), ..., (%g,%g)\n",
 - 
 
898 
			xxx[0],yyy[0], xxx[1],yyy[1], xxx[ns-1],yyy[ns-1]);
 - 
 
899 
	    else
 - 
 
900 
		for(iii = 0; iii < ns; iii++)
 - 
 
901 
		    Rprintf(" (%g,%g)%s", xxx[iii],yyy[iii],
 - 
 
902 
			    (iii < ns-1) ? "," : "\n");
 - 
 
903 
#endif
 - 
 
904 
 
 - 
 
905 
//	    GMode(1, dd);
 - 
 
906 
 
 - 
 
907 
	    if (drawLabels) {
 - 
 
908 
		/* If user supplied labels, use i'th one of them
 - 
 
909 
		   Otherwise stringify the z-value of the contour */
 - 
 
910 
		cetype_t enc = CE_NATIVE;
 - 
 
911 
		buffer[0] = ' ';
 - 
 
912 
		if (!isNull(labels)) {
 - 
 
913 
		    int numl = length(labels);
 - 
 
914 
		    strcpy(&buffer[1], CHAR(STRING_ELT(labels, cnum % numl)));
 - 
 
915 
		    enc = getCharCE(STRING_ELT(labels, cnum % numl));
 - 
 
916 
		}
 - 
 
917 
		else {
 - 
 
918 
		    PROTECT(lab = allocVector(REALSXP, 1));
 - 
 
919 
		    REAL(lab)[0] = zc;
 - 
 
920 
		    lab = labelformat(lab);
 - 
 
921 
		    strcpy(&buffer[1], CHAR(STRING_ELT(lab, 0))); /* ASCII */
 - 
 
922 
		    UNPROTECT(1);
 - 
 
923 
		}
 - 
 
924 
		buffer[strlen(buffer)+1] = '\0';
 - 
 
925 
		buffer[strlen(buffer)] = ' ';
 - 
 
926 
 
 - 
 
927 
		labelDistance = GStrWidth(buffer, enc, INCHES, dd);
 - 
 
928 
		labelHeight = GStrHeight(buffer, enc, INCHES, dd);
 - 
 
929 
 
 - 
 
930 
		if (labelDistance > 0) {
 - 
 
931 
		    /* Try to find somewhere to draw the label */
 - 
 
932 
		    switch (method) {
 - 
 
933 
		    case 0: /* draw label at one end of contour
 - 
 
934 
			       overwriting contour line
 - 
 
935 
			    */
 - 
 
936 
			if (useStart(xxx, yyy, ns, dd) )
 - 
 
937 
			    indx = 0;
 - 
 
938 
			else
 - 
 
939 
			    indx = ns - 1;
 - 
 
940 
			break;
 - 
 
941 
		    case 1: /* draw label at one end of contour
 - 
 
942 
			       embedded in contour
 - 
 
943 
			       no overlapping labels
 - 
 
944 
			    */
 - 
 
945 
			indx = 0;
 - 
 
946 
			range = 0;
 - 
 
947 
			gotLabel = FALSE;
 - 
 
948 
			if (useStart(xxx, yyy, ns, dd)) {
 - 
 
949 
			    iii = 0;
 - 
 
950 
			    n = findGapUp(xxx, yyy, ns, labelDistance, dd);
 - 
 
951 
			}
 - 
 
952 
			else {
 - 
 
953 
			    n = findGapDown(xxx, yyy, ns, labelDistance, dd);
 - 
 
954 
			    iii = ns - n - 1;
 - 
 
955 
			}
 - 
 
956 
			if (n > 0) {
 - 
 
957 
			    /** Find 4 corners of label extents **/
 - 
 
958 
			    FindCorners(labelDistance, labelHeight, label1,
 - 
 
959 
					xxx[iii], yyy[iii],
 - 
 
960 
					xxx[iii+n], yyy[iii+n], dd);
 - 
 
961 
 
 - 
 
962 
			    /** Test corners for intersection with previous labels **/
 - 
 
963 
			    label2 = labelList;
 - 
 
964 
			    result = 0;
 - 
 
965 
			    while ((result == 0) && (label2 != R_NilValue)) {
 - 
 
966 
				result = TestLabelIntersection(label1, CAR(label2));
 - 
 
967 
				label2 = CDR(label2);
 - 
 
968 
			    }
 - 
 
969 
			    if (result == 0) {
 - 
 
970 
				result = LabelInsideWindow(label1, dd);
 - 
 
971 
				if (result == 0) {
 - 
 
972 
				    indx = iii;
 - 
 
973 
				    range = n;
 - 
 
974 
				    gotLabel = TRUE;
 - 
 
975 
				}
 - 
 
976 
			    }
 - 
 
977 
			}
 - 
 
978 
			break;
 - 
 
979 
		    case 2: /* draw label on flattest portion of contour
 - 
 
980 
			       embedded in contour line
 - 
 
981 
			       no overlapping labels
 - 
 
982 
			    */
 - 
 
983 
			/* Look for flatest sequence of contour gradients */
 - 
 
984 
			lowestVariance = 9999999;   /* A large number */
 - 
 
985 
			indx = 0;
 - 
 
986 
			range = 0;
 - 
 
987 
			gotLabel = FALSE;
 - 
 
988 
			for (iii = 0; iii < ns; iii++) {
 - 
 
989 
			    distanceSum = 0;
 - 
 
990 
			    avgGradient = 0;
 - 
 
991 
			    squareSum = 0;
 - 
 
992 
			    n = 0;
 - 
 
993 
			    jjj = (iii + 1);
 - 
 
994 
			    while ((jjj < ns-1) &&
 - 
 
995 
				   (distanceSum < labelDistance)) {
 - 
 
996 
 
 - 
 
997 
				/* Find a gap big enough for the label
 - 
 
998 
				   use several segments if necessary
 - 
 
999 
				*/
 - 
 
1000 
				dX = xxx[jjj] - xxx[jjj - n - 1];
 - 
 
1001 
				dY = yyy[jjj] - yyy[jjj - n - 1];
 - 
 
1002 
				dXC = GConvertXUnits(dX, USER, INCHES, dd);
 - 
 
1003 
				dYC = GConvertYUnits(dY, USER, INCHES, dd);
 - 
 
1004 
				distanceSum = hypot(dXC, dYC);
 - 
 
1005 
 
 - 
 
1006 
				/* Calculate the variance of the gradients
 - 
 
1007 
				   of the segments that will make way for the
 - 
 
1008 
				   label
 - 
 
1009 
				*/
 - 
 
1010 
				deltaX = xxx[jjj] - xxx[jjj - 1];
 - 
 
1011 
				deltaY = yyy[jjj] - yyy[jjj - 1];
 - 
 
1012 
				if (deltaX == 0) {deltaX = 1;}
 - 
 
1013 
				avgGradient += (deltaY/deltaX);
 - 
 
1014 
				squareSum += avgGradient * avgGradient;
 - 
 
1015 
				jjj = (jjj + 1);
 - 
 
1016 
				n += 1;
 - 
 
1017 
			    }
 - 
 
1018 
			    if (distanceSum < labelDistance)
 - 
 
1019 
				break;
 - 
 
1020 
 
 - 
 
1021 
			    /** Find 4 corners of label extents **/
 - 
 
1022 
			    FindCorners(labelDistance, labelHeight, label1,
 - 
 
1023 
					xxx[iii], yyy[iii],
 - 
 
1024 
					xxx[iii+n], yyy[iii+n], dd);
 - 
 
1025 
 
 - 
 
1026 
			    /** Test corners for intersection with previous labels **/
 - 
 
1027 
			    label2 = labelList;
 - 
 
1028 
			    result = 0;
 - 
 
1029 
			    while ((result == 0) && (label2 != R_NilValue)) {
 - 
 
1030 
				result = TestLabelIntersection(label1, CAR(label2));
 - 
 
1031 
				label2 = CDR(label2);
 - 
 
1032 
			    }
 - 
 
1033 
			    if (result == 0)
 - 
 
1034 
				result = LabelInsideWindow(label1, dd);
 - 
 
1035 
			    if (result == 0) {
 - 
 
1036 
				variance = (squareSum - (avgGradient * avgGradient) / n) / n;
 - 
 
1037 
				avgGradient /= n;
 - 
 
1038 
				if (variance < lowestVariance) {
 - 
 
1039 
				    lowestVariance = variance;
 - 
 
1040 
				    indx = iii;
 - 
 
1041 
				    range = n;
 - 
 
1042 
				}
 - 
 
1043 
			    }
 - 
 
1044 
			    if (lowestVariance < 9999999)
 - 
 
1045 
				gotLabel = TRUE;
 - 
 
1046 
			}
 - 
 
1047 
		    } /* switch (method) */
 - 
 
1048 
 
 - 
 
1049 
		    if (method == 0) {
 - 
 
1050 
			GPolyline(ns, xxx, yyy, USER, dd);
 - 
 
1051 
			GText(xxx[indx], yyy[indx], USER, buffer,
 - 
 
1052 
			      CE_NATIVE/*FIX*/,
 - 
 
1053 
			      .5, .5, 0, dd);
 - 
 
1054 
		    }
 - 
 
1055 
		    else {
 - 
 
1056 
			if (indx > 0)
 - 
 
1057 
		            GPolyline(indx+1, xxx, yyy, USER, dd);
 - 
 
1058 
			if (ns-1-indx-range > 0)
 - 
 
1059 
			    GPolyline(ns-indx-range, xxx+indx+range, yyy+indx+range,
 - 
 
1060 
			          USER, dd);
 - 
 
1061 
			if (gotLabel) {
 - 
 
1062 
			    /* find which plot edge we are closest to */
 - 
 
1063 
			    int closest; /* 0 = indx,  1 = indx+range */
 - 
 
1064 
			    double dx1, dx2, dy1, dy2, dmin;
 - 
 
1065 
			    dx1 = fmin2((xxx[indx] - gpptr(dd)->usr[0]),
 - 
 
1066 
					(gpptr(dd)->usr[1] - xxx[indx]));
 - 
 
1067 
			    dx2 = fmin2((gpptr(dd)->usr[1] - xxx[indx+range]),
 - 
 
1068 
					(xxx[indx+range] - gpptr(dd)->usr[0]));
 - 
 
1069 
			    if (dx1 < dx2) {
 - 
 
1070 
				closest = 0;
 - 
 
1071 
				dmin = dx1;
 - 
 
1072 
			    } else {
 - 
 
1073 
				closest = 1;
 - 
 
1074 
				dmin = dx2;
 - 
 
1075 
			    }
 - 
 
1076 
			    dy1 = fmin2((yyy[indx] - gpptr(dd)->usr[2]),
 - 
 
1077 
					(gpptr(dd)->usr[3] - yyy[indx]));
 - 
 
1078 
			    if (closest && (dy1 < dmin)) {
 - 
 
1079 
				closest = 0;
 - 
 
1080 
				dmin = dy1;
 - 
 
1081 
			    } else if (dy1 < dmin)
 - 
 
1082 
				dmin = dy1;
 - 
 
1083 
			    dy2 = fmin2((gpptr(dd)->usr[3] - yyy[indx+range]),
 - 
 
1084 
					(yyy[indx+range] - gpptr(dd)->usr[2]));
 - 
 
1085 
			    if (!closest && (dy2 < dmin))
 - 
 
1086 
				closest = 1;
 - 
 
1087 
 
 - 
 
1088 
			    dx = GConvertXUnits(xxx[indx+range] - xxx[indx],
 - 
 
1089 
						USER, INCHES, dd);
 - 
 
1090 
			    dy = GConvertYUnits(yyy[indx+range] - yyy[indx],
 - 
 
1091 
						USER, INCHES, dd);
 - 
 
1092 
			    dxy = hypot(dx, dy);
 - 
 
1093 
 
 - 
 
1094 
			    /* save the current label for checking overlap */
 - 
 
1095 
			    label2 = allocVector(REALSXP, 8);
 - 
 
1096 
 
 - 
 
1097 
			    FindCorners(labelDistance, labelHeight, label2,
 - 
 
1098 
					xxx[indx], yyy[indx],
 - 
 
1099 
					xxx[indx+range], yyy[indx+range], dd);
 - 
 
1100 
			    UNPROTECT_PTR(labelList);
 - 
 
1101 
			    labelList = PROTECT(CONS(label2, labelList));
 - 
 
1102 
 
 - 
 
1103 
			    ddl = FALSE;
 - 
 
1104 
			    /* draw an extra bit of segment if the label
 - 
 
1105 
			       doesn't fill the gap */
 - 
 
1106 
			    if (closest) {
 - 
 
1107 
				xStart = xxx[indx+range] -
 - 
 
1108 
				    (xxx[indx+range] - xxx[indx]) *
 - 
 
1109 
				    labelDistance / dxy;
 - 
 
1110 
				yStart = yyy[indx+range] -
 - 
 
1111 
				    (yyy[indx+range] - yyy[indx]) *
 - 
 
1112 
				    labelDistance / dxy;
 - 
 
1113 
				if (labelDistance / dxy < 1)
 - 
 
1114 
				    GLine(xxx[indx], yyy[indx],
 - 
 
1115 
					  xStart, yStart,
 - 
 
1116 
					  USER, dd);
 - 
 
1117 
			    } else {
 - 
 
1118 
				xStart = xxx[indx] +
 - 
 
1119 
				    (xxx[indx+range] - xxx[indx]) *
 - 
 
1120 
				    labelDistance / dxy;
 - 
 
1121 
				yStart = yyy[indx] +
 - 
 
1122 
				    (yyy[indx+range] - yyy[indx]) *
 - 
 
1123 
				    labelDistance / dxy;
 - 
 
1124 
				if (labelDistance / dxy < 1)
 - 
 
1125 
				    GLine(xStart, yStart,
 - 
 
1126 
					  xxx[indx+range], yyy[indx+range],
 - 
 
1127 
					  USER, dd);
 - 
 
1128 
			    }
 - 
 
1129 
 
 - 
 
1130 
			    /*** Draw contour labels ***/
 - 
 
1131 
			    if (xxx[indx] < xxx[indx+range]) {
 - 
 
1132 
				if (closest) {
 - 
 
1133 
				    ux = xStart;
 - 
 
1134 
				    uy = yStart;
 - 
 
1135 
				    vx = xxx[indx+range];
 - 
 
1136 
				    vy = yyy[indx+range];
 - 
 
1137 
				} else {
 - 
 
1138 
				    ux = xxx[indx];
 - 
 
1139 
				    uy = yyy[indx];
 - 
 
1140 
				    vx = xStart;
 - 
 
1141 
				    vy = yStart;
 - 
 
1142 
				}
 - 
 
1143 
			    }
 - 
 
1144 
			    else {
 - 
 
1145 
				if (closest) {
 - 
 
1146 
				    ux = xxx[indx+range];
 - 
 
1147 
				    uy = yyy[indx+range];
 - 
 
1148 
				    vx = xStart;
 - 
 
1149 
				    vy = yStart;
 - 
 
1150 
				} else {
 - 
 
1151 
				    ux = xStart;
 - 
 
1152 
				    uy = yStart;
 - 
 
1153 
				    vx = xxx[indx];
 - 
 
1154 
				    vy = yyy[indx];
 - 
 
1155 
				}
 - 
 
1156 
			    }
 - 
 
1157 
 
 - 
 
1158 
			    if (!ddl) {
 - 
 
1159 
				/* convert to INCHES for calculation of
 - 
 
1160 
				   angle to draw text
 - 
 
1161 
				*/
 - 
 
1162 
				GConvert(&ux, &uy, USER, INCHES, dd);
 - 
 
1163 
				GConvert(&vx, &vy, USER, INCHES, dd);
 - 
 
1164 
				/* 0, .5 => left, centre justified */
 - 
 
1165 
				GText (ux, uy, INCHES, buffer,
 - 
 
1166 
				       CE_NATIVE/*FIX*/,0, .5,
 - 
 
1167 
				       (180 / 3.14) * atan2(vy - uy, vx - ux),
 - 
 
1168 
				       dd);
 - 
 
1169 
			    }
 - 
 
1170 
			} /* if (gotLabel) */
 - 
 
1171 
		    } /* if (method == 0) else ... */
 - 
 
1172 
		} /* if (labelDistance > 0) */
 - 
 
1173 
 
 - 
 
1174 
	    } /* if (drawLabels) */
 - 
 
1175 
	    else {
 - 
 
1176 
		GPolyline(ns, xxx, yyy, USER, dd);
 - 
 
1177 
	    }
 - 
 
1178 
 
 - 
 
1179 
//	    GMode(0, dd);
 - 
 
1180 
	    vmaxset(vmax);
 - 
 
1181 
	} /* while */
 - 
 
1182 
      } /* for(i .. )  for(j ..) */
 - 
 
1183 
    vmaxset(vmax); /* now we are done with ctr_SegDB */
 - 
 
1184 
    UNPROTECT_PTR(label1); /* pwwwargh! This is messy, but last thing
 - 
 
1185 
			      protected is likely labelList, and that needs
 - 
 
1186 
			      to be preserved across calls */
 - 
 
1187 
}
 - 
 
1188 
 
 - 
 
1189 
 
 - 
 
1190 
SEXP attribute_hidden Rg_contourDef(void)
 - 
 
1191 
{
 - 
 
1192 
    return ScalarLogical(GEcurrentDevice()->dev->useRotatedTextInContour);
 - 
 
1193 
}
 - 
 
1194 
 
 - 
 
1195 
/* contour(x, y, z, levels, labels, labcex, drawlabels,
 - 
 
1196 
 *         method, vfont, col = col, lty = lty, lwd = lwd)
 - 
 
1197 
 */
 - 
 
1198 
SEXP attribute_hidden do_contour(SEXP call, SEXP op, SEXP args, SEXP env)
 - 
 
1199 
{
 - 
 
1200 
    SEXP oargs, c, x, y, z, vfont, col, rawcol, lty, lwd, labels;
 - 
 
1201 
    int i, j, nx, ny, nc, ncol, nlty, nlwd;
 - 
 
1202 
    int ltysave, fontsave = 1 /* -Wall */;
 - 
 
1203 
    rcolor colsave;
 - 
 
1204 
    double cexsave, lwdsave;
 - 
 
1205 
    double atom, zmin, zmax;
 - 
 
1206 
    const void *vmax, *vmax0;
 - 
 
1207 
    char familysave[201];
 - 
 
1208 
    int method;
 - 
 
1209 
    Rboolean drawLabels;
 - 
 
1210 
    double labcex;
 - 
 
1211 
    pGEDevDesc dd = GEcurrentDevice();
 - 
 
1212 
    SEXP result = R_NilValue;
 - 
 
1213 
 
 - 
 
1214 
    GCheckState(dd);
 - 
 
1215 
 
 - 
 
1216 
    if (length(args) < 4)
 - 
 
1217 
	error(_("too few arguments"));
 - 
 
1218 
    PrintDefaults(); /* prepare for labelformat */
 - 
 
1219 
 
 - 
 
1220 
    oargs = args;
 - 
 
1221 
 
 - 
 
1222 
    x = CAR(args);
 - 
 
1223 
    internalTypeCheck(call, x, REALSXP);
 - 
 
1224 
    nx = LENGTH(x);
 - 
 
1225 
    args = CDR(args);
 - 
 
1226 
 
 - 
 
1227 
    y = CAR(args);
 - 
 
1228 
    internalTypeCheck(call, y, REALSXP);
 - 
 
1229 
    ny = LENGTH(y);
 - 
 
1230 
    args = CDR(args);
 - 
 
1231 
 
 - 
 
1232 
    z = CAR(args);
 - 
 
1233 
    internalTypeCheck(call, z, REALSXP);
 - 
 
1234 
    args = CDR(args);
 - 
 
1235 
 
 - 
 
1236 
    /* levels */
 - 
 
1237 
    c = CAR(args);
 - 
 
1238 
    internalTypeCheck(call, c, REALSXP);
 - 
 
1239 
    nc = LENGTH(c);
 - 
 
1240 
    args = CDR(args);
 - 
 
1241 
 
 - 
 
1242 
    labels = CAR(args);
 - 
 
1243 
    if (!isNull(labels))
 - 
 
1244 
	internalTypeCheck(call, labels, STRSXP);
 - 
 
1245 
    args = CDR(args);
 - 
 
1246 
 
 - 
 
1247 
    labcex = asReal(CAR(args));
 - 
 
1248 
    args = CDR(args);
 - 
 
1249 
 
 - 
 
1250 
    drawLabels = (Rboolean)asLogical(CAR(args));
 - 
 
1251 
    args = CDR(args);
 - 
 
1252 
 
 - 
 
1253 
    method = asInteger(CAR(args)); args = CDR(args);
 - 
 
1254 
    if (method < 1 || method > 3)
 - 
 
1255 
	error(_("invalid '%s' value"), "method");
 - 
 
1256 
 
 - 
 
1257 
    PROTECT(vfont = FixupVFont(CAR(args)));
 - 
 
1258 
    if (!isNull(vfont)) {
 - 
 
1259 
	strncpy(familysave, gpptr(dd)->family, 201);
 - 
 
1260 
	strncpy(gpptr(dd)->family, "Her ", 201);
 - 
 
1261 
	gpptr(dd)->family[3] = (char) INTEGER(vfont)[0];
 - 
 
1262 
	fontsave = gpptr(dd)->font;
 - 
 
1263 
	gpptr(dd)->font = INTEGER(vfont)[1];
 - 
 
1264 
    }
 - 
 
1265 
    args = CDR(args);
 - 
 
1266 
 
 - 
 
1267 
    rawcol = CAR(args);
 - 
 
1268 
    PROTECT(col = FixupCol(rawcol, R_TRANWHITE));
 - 
 
1269 
    ncol = length(col);
 - 
 
1270 
    args = CDR(args);
 - 
 
1271 
 
 - 
 
1272 
    PROTECT(lty = FixupLty(CAR(args), gpptr(dd)->lty));
 - 
 
1273 
    nlty = length(lty);
 - 
 
1274 
    args = CDR(args);
 - 
 
1275 
 
 - 
 
1276 
    PROTECT(lwd = FixupLwd(CAR(args), gpptr(dd)->lwd));
 - 
 
1277 
    nlwd = length(lwd);
 - 
 
1278 
    args = CDR(args);
 - 
 
1279 
 
 - 
 
1280 
    if (nx < 2 || ny < 2)
 - 
 
1281 
	error(_("insufficient 'x' or 'y' values"));
 - 
 
1282 
 
 - 
 
1283 
    if (nrows(z) != nx || ncols(z) != ny)
 - 
 
1284 
	error(_("dimension mismatch"));
 - 
 
1285 
 
 - 
 
1286 
    if (nc < 1)
 - 
 
1287 
	error(_("no contour values"));
 - 
 
1288 
 
 - 
 
1289 
    for (i = 0; i < nx; i++) {
 - 
 
1290 
	if (!R_FINITE(REAL(x)[i]))
 - 
 
1291 
	    error(_("missing 'x' values"));
 - 
 
1292 
	if (i > 0 && REAL(x)[i] < REAL(x)[i - 1])
 - 
 
1293 
	    error(_("increasing 'x' values expected"));
 - 
 
1294 
    }
 - 
 
1295 
 
 - 
 
1296 
    for (i = 0; i < ny; i++) {
 - 
 
1297 
	if (!R_FINITE(REAL(y)[i]))
 - 
 
1298 
	    error(_("missing 'y' values"));
 - 
 
1299 
	if (i > 0 && REAL(y)[i] < REAL(y)[i - 1])
 - 
 
1300 
	    error(_("increasing 'y' values expected"));
 - 
 
1301 
    }
 - 
 
1302 
 
 - 
 
1303 
    for (i = 0; i < nc; i++)
 - 
 
1304 
	if (!R_FINITE(REAL(c)[i]))
 - 
 
1305 
	    error(_("invalid NA contour values"));
 - 
 
1306 
 
 - 
 
1307 
    zmin = DBL_MAX;
 - 
 
1308 
    zmax = DBL_MIN;
 - 
 
1309 
    for (i = 0; i < nx * ny; i++)
 - 
 
1310 
	if (R_FINITE(REAL(z)[i])) {
 - 
 
1311 
	    if (zmax < REAL(z)[i]) zmax =  REAL(z)[i];
 - 
 
1312 
	    if (zmin > REAL(z)[i]) zmin =  REAL(z)[i];
 - 
 
1313 
	}
 - 
 
1314 
 
 - 
 
1315 
    if (zmin >= zmax) {
 - 
 
1316 
	if (zmin == zmax)
 - 
 
1317 
	    warning(_("all z values are equal"));
 - 
 
1318 
	else
 - 
 
1319 
	    warning(_("all z values are NA"));
 - 
 
1320 
	UNPROTECT(4);
 - 
 
1321 
	return R_NilValue;
 - 
 
1322 
    }
 - 
 
1323 
 
 - 
 
1324 
    /* change to 1e-3, reconsidered because of PR#897
 - 
 
1325 
     * but 1e-7, and even  2*DBL_EPSILON do not prevent inf.loop in contour().
 - 
 
1326 
     * maybe something like   16 * DBL_EPSILON * (..).
 - 
 
1327 
     * see also max_contour_segments above */
 - 
 
1328 
    atom = 1e-3 * (zmax - zmin);
 - 
 
1329 
 
 - 
 
1330 
    /* Initialize the segment data base */
 - 
 
1331 
 
 - 
 
1332 
    /* Note we must be careful about resetting */
 - 
 
1333 
    /* the top of the stack, otherwise we run out of */
 - 
 
1334 
    /* memory after a sequence of displaylist replays */
 - 
 
1335 
 
 - 
 
1336 
    vmax0 = vmaxget();
 - 
 
1337 
    ctr_SegDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));
 - 
 
1338 
 
 - 
 
1339 
    for (i = 0; i < nx; i++)
 - 
 
1340 
	for (j = 0; j < ny; j++)
 - 
 
1341 
	    ctr_SegDB[i + j * nx] = NULL;
 - 
 
1342 
 
 - 
 
1343 
    /* Draw the contours -- note the heap release */
 - 
 
1344 
 
 - 
 
1345 
    ltysave = gpptr(dd)->lty;
 - 
 
1346 
    colsave = gpptr(dd)->col;
 - 
 
1347 
    lwdsave = gpptr(dd)->lwd;
 - 
 
1348 
    cexsave = gpptr(dd)->cex;
 - 
 
1349 
    labelList = PROTECT(R_NilValue);
 - 
 
1350 
 
 - 
 
1351 
 
 - 
 
1352 
    /* draw contour for levels[i] */
 - 
 
1353 
    GMode(1, dd);
 - 
 
1354 
    for (i = 0; i < nc; i++) {
 - 
 
1355 
	vmax = vmaxget();
 - 
 
1356 
	gpptr(dd)->lty = INTEGER(lty)[i % nlty];
 - 
 
1357 
	if (gpptr(dd)->lty == NA_INTEGER)
 - 
 
1358 
	    gpptr(dd)->lty = ltysave;
 - 
 
1359 
	if (isNAcol(rawcol, i, ncol))
 - 
 
1360 
	    gpptr(dd)->col = colsave;
 - 
 
1361 
	else
 - 
 
1362 
	    gpptr(dd)->col = INTEGER(col)[i % ncol];
 - 
 
1363 
	gpptr(dd)->lwd = REAL(lwd)[i % nlwd];
 - 
 
1364 
	if (!R_FINITE(gpptr(dd)->lwd))
 - 
 
1365 
	    gpptr(dd)->lwd = lwdsave;
 - 
 
1366 
	gpptr(dd)->cex = labcex;
 - 
 
1367 
	contour(x, nx, y, ny, z, REAL(c)[i], labels, i,
 - 
 
1368 
		drawLabels, method - 1, atom, dd);
 - 
 
1369 
	vmaxset(vmax);
 - 
 
1370 
    }
 - 
 
1371 
    GMode(0, dd);
 - 
 
1372 
    vmaxset(vmax0);
 - 
 
1373 
    gpptr(dd)->lty = ltysave;
 - 
 
1374 
    gpptr(dd)->col = colsave;
 - 
 
1375 
    gpptr(dd)->lwd = lwdsave;
 - 
 
1376 
    gpptr(dd)->cex = cexsave;
 - 
 
1377 
    if(!isNull(vfont)) {
 - 
 
1378 
	strncpy(gpptr(dd)->family, familysave, 201);
 - 
 
1379 
	gpptr(dd)->font = fontsave;
 - 
 
1380 
    }
 - 
 
1381 
    UNPROTECT(5);
 - 
 
1382 
    /* NOTE: only record operation if no "error"  */
 - 
 
1383 
    /* NOTE: on replay, call == R_NilValue */
 - 
 
1384 
    if (GRecording(call, dd))
 - 
 
1385 
	GErecordGraphicOperation(op, oargs, dd);
 - 
 
1386 
    return result;
 - 
 
1387 
}
 -