1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
4 * Copyright (C) 1997--2013 The R Core Team
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, a copy is available at
18 * https://www.R-project.org/Licenses/
19 */
20
21
22 /* Stuff for labels on contour plots
23 Originally written by Nicholas Hildreth
24 Adapted by Paul Murrell
25 */
26
27 /* Included by src/main/plot3d.c and src/library/graphics/src/plot3d */
28
29 /* C o n t o u r P l o t t i n g */
30
31 typedef struct SEG {
32 struct SEG *next;
33 double x0;
34 double y0;
35 double x1;
36 double y1;
37 } SEG, *SEGP;
38
39
ctr_intersect(double z0,double z1,double zc,double * f)40 static int ctr_intersect(double z0, double z1, double zc, double *f)
41 {
42 /* Old test was ((z0 - zc) * (z1 - zc) < 0.0), but rounding led to inconsistencies
43 in PR#15454 */
44 if ( (z0 < zc) != (z1 < zc) && z0 != zc && z1 != zc ) {
45 *f = (zc - z0) / (z1 - z0);
46 return 1;
47 }
48 return 0;
49 }
50
ctr_newseg(double x0,double y0,double x1,double y1,SEGP prev)51 static SEGP ctr_newseg(double x0, double y0, double x1, double y1, SEGP prev)
52 {
53 SEGP seg = (SEGP)R_alloc(1, sizeof(SEG));
54 seg->x0 = x0;
55 seg->y0 = y0;
56 seg->x1 = x1;
57 seg->y1 = y1;
58 seg->next = prev;
59 return seg;
60 }
61
ctr_swapseg(SEGP seg)62 static void ctr_swapseg(SEGP seg)
63 {
64 double x, y;
65 x = seg->x0;
66 y = seg->y0;
67 seg->x0 = seg->x1;
68 seg->y0 = seg->y1;
69 seg->x1 = x;
70 seg->y1 = y;
71 }
72
73 /* ctr_segdir(): Determine the entry direction to the next cell */
74 /* and update the cell indices */
75
76 #define XMATCH(x0,x1) (fabs(x0-x1) == 0)
77 #define YMATCH(y0,y1) (fabs(y0-y1) == 0)
78
ctr_segdir(double xend,double yend,double * x,double * y,int * i,int * j,int nx,int ny)79 static int ctr_segdir(double xend, double yend, double *x, double *y,
80 int *i, int *j, int nx, int ny)
81 {
82 if (YMATCH(yend, y[*j])) {
83 if (*j == 0)
84 return 0;
85 *j = *j - 1;
86 return 3;
87 }
88 if (XMATCH(xend, x[*i])) {
89 if (*i == 0)
90 return 0;
91 *i = *i - 1;
92 return 4;
93 }
94 if (YMATCH(yend, y[*j + 1])) {
95 if (*j >= ny - 1)
96 return 0;
97 *j = *j + 1;
98 return 1;
99 }
100 if (XMATCH(xend, x[*i + 1])) {
101 if (*i >= nx - 1)
102 return 0;
103 *i = *i + 1;
104 return 2;
105 }
106 return 0;
107 }
108
109 /* Search seglist for a segment with endpoint (xend, yend). */
110 /* The cell entry direction is dir, and if tail=1/0 we are */
111 /* building the tail/head of a contour. The matching segment */
112 /* is pointed to by seg and the updated segment list (with */
113 /* the matched segment stripped) is returned by the funtion. */
114
ctr_segupdate(double xend,double yend,int dir,Rboolean tail,SEGP seglist,SEGP * seg)115 static SEGP ctr_segupdate(double xend, double yend, int dir, Rboolean tail,
116 SEGP seglist, SEGP* seg)
117 {
118 if (seglist == NULL) {
119 *seg = NULL;
120 return NULL;
121 }
122 switch (dir) {
123 case 1:
124 case 3:
125 if (YMATCH(yend,seglist->y0)) {
126 if (!tail)
127 ctr_swapseg(seglist);
128 *seg = seglist;
129 return seglist->next;
130 }
131 if (YMATCH(yend,seglist->y1)) {
132 if (tail)
133 ctr_swapseg(seglist);
134 *seg = seglist;
135 return seglist->next;
136 }
137 break;
138 case 2:
139 case 4:
140 if (XMATCH(xend,seglist->x0)) {
141 if (!tail)
142 ctr_swapseg(seglist);
143 *seg = seglist;
144 return seglist->next;
145 }
146 if (XMATCH(xend,seglist->x1)) {
147 if (tail)
148 ctr_swapseg(seglist);
149 *seg = seglist;
150 return seglist->next;
151 }
152 break;
153 }
154 seglist->next = ctr_segupdate(xend, yend, dir, tail, seglist->next, seg);
155 return seglist;
156 }
157
158
159
160 /*
161 * Generate a list of segments for a single level
162 *
163 * NB this R_allocs its return value, so callers need to manage R_alloc stack.
164 */
contourLines(double * x,int nx,double * y,int ny,double * z,double zc,double atom)165 static SEGP* contourLines(double *x, int nx, double *y, int ny,
166 double *z, double zc, double atom)
167 {
168 double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];
169 int i, j, k, l, m, nacode;
170 SEGP seglist;
171 SEGP *segmentDB;
172 /* Initialize the segment data base */
173 /* Note we must be careful about resetting */
174 /* the top of the stack, otherwise we run out of */
175 /* memory after a sequence of displaylist replays */
176 /*
177 * This reset is done out in GEcontourLines
178 */
179 segmentDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));
180 for (i = 0; i < nx; i++)
181 for (j = 0; j < ny; j++)
182 segmentDB[i + j * nx] = NULL;
183 for (i = 0; i < nx - 1; i++) {
184 xl = x[i];
185 xh = x[i + 1];
186 for (j = 0; j < ny - 1; j++) {
187 yl = y[j];
188 yh = y[j + 1];
189 k = i + j * nx;
190 zll = z[k];
191 zhl = z[k + 1];
192 zlh = z[k + nx];
193 zhh = z[k + nx + 1];
194
195 /* If the value at a corner is exactly equal to a contour level,
196 * change that value by a tiny amount */
197
198 if (zll == zc) zll += atom;
199 if (zhl == zc) zhl += atom;
200 if (zlh == zc) zlh += atom;
201 if (zhh == zc) zhh += atom;
202 #ifdef DEBUG_contour
203 /* Haven't seen this happening (MM): */
204 if (zll == zc) REprintf(" [%d,%d] ll: %g\n",i,j, zll);
205 if (zhl == zc) REprintf(" [%d,%d] hl: %g\n",i,j, zhl);
206 if (zlh == zc) REprintf(" [%d,%d] lh: %g\n",i,j, zlh);
207 if (zhh == zc) REprintf(" [%d,%d] hh: %g\n",i,j, zhh);
208 #endif
209 /* Check for intersections with sides */
210
211 nacode = 0;
212 if (R_FINITE(zll)) nacode += 1;
213 if (R_FINITE(zhl)) nacode += 2;
214 if (R_FINITE(zlh)) nacode += 4;
215 if (R_FINITE(zhh)) nacode += 8;
216
217 k = 0;
218 switch (nacode) {
219 case 15:
220 if (ctr_intersect(zll, zhl, zc, &f)) {
221 xx[k] = xl + f * (xh - xl);
222 yy[k] = yl; k++;
223 }
224 if (ctr_intersect(zll, zlh, zc, &f)) {
225 yy[k] = yl + f * (yh - yl);
226 xx[k] = xl; k++;
227 }
228 if (ctr_intersect(zhl, zhh, zc, &f)) {
229 yy[k] = yl + f * (yh - yl);
230 xx[k] = xh; k++;
231 }
232 if (ctr_intersect(zlh, zhh, zc, &f)) {
233 xx[k] = xl + f * (xh - xl);
234 yy[k] = yh; k++;
235 }
236 break;
237 case 14:
238 if (ctr_intersect(zhl, zhh, zc, &f)) {
239 yy[k] = yl + f * (yh - yl);
240 xx[k] = xh; k++;
241 }
242 if (ctr_intersect(zlh, zhh, zc, &f)) {
243 xx[k] = xl + f * (xh - xl);
244 yy[k] = yh; k++;
245 }
246 if (ctr_intersect(zlh, zhl, zc, &f)) {
247 xx[k] = xl + f * (xh - xl);
248 yy[k] = yh + f * (yl - yh);
249 k++;
250 }
251 break;
252 case 13:
253 if (ctr_intersect(zll, zlh, zc, &f)) {
254 yy[k] = yl + f * (yh - yl);
255 xx[k] = xl; k++;
256 }
257 if (ctr_intersect(zlh, zhh, zc, &f)) {
258 xx[k] = xl + f * (xh - xl);
259 yy[k] = yh; k++;
260 }
261 if (ctr_intersect(zll, zhh, zc, &f)) {
262 xx[k] = xl + f * (xh - xl);
263 yy[k] = yl + f * (yh - yl);
264 k++;
265 }
266 break;
267 case 11:
268 if (ctr_intersect(zhl, zhh, zc, &f)) {
269 yy[k] = yl + f * (yh - yl);
270 xx[k] = xh; k++;
271 }
272 if (ctr_intersect(zll, zhl, zc, &f)) {
273 xx[k] = xl + f * (xh - xl);
274 yy[k] = yl; k++;
275 }
276 if (ctr_intersect(zll, zhh, zc, &f)) {
277 xx[k] = xl + f * (xh - xl);
278 yy[k] = yl + f * (yh - yl);
279 k++;
280 }
281 break;
282 case 7:
283 if (ctr_intersect(zll, zlh, zc, &f)) {
284 yy[k] = yl + f * (yh - yl);
285 xx[k] = xl; k++;
286 }
287 if (ctr_intersect(zll, zhl, zc, &f)) {
288 xx[k] = xl + f * (xh - xl);
289 yy[k] = yl; k++;
290 }
291 if (ctr_intersect(zlh, zhl, zc, &f)) {
292 xx[k] = xl + f * (xh - xl);
293 yy[k] = yh + f * (yl - yh);
294 k++;
295 }
296 break;
297 }
298
299 /* We now have k(=2,4) endpoints */
300 /* Decide which to join */
301
302 seglist = NULL;
303
304 if (k > 0) {
305 if (k == 2) {
306 seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
307 }
308 else if (k == 4) {
309 for (k = 3; k >= 1; k--) {
310 m = k;
311 xl = xx[k];
312 for (l = 0; l < k; l++) {
313 if (xx[l] > xl) {
314 xl = xx[l];
315 m = l;
316 }
317 }
318 if (m != k) {
319 xl = xx[k];
320 yl = yy[k];
321 xx[k] = xx[m];
322 yy[k] = yy[m];
323 xx[m] = xl;
324 yy[m] = yl;
325 }
326 }
327 seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
328 seglist = ctr_newseg(xx[2], yy[2], xx[3], yy[3], seglist);
329 }
330 else error("k = %d, should be 2 or 4", k);
331 }
332 segmentDB[i + j * nx] = seglist;
333 }
334 }
335 return segmentDB;
336 }
337