1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
4 * Copyright (C) 1997--2014 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 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 #include <Defn.h>
27 #include <Internal.h>
28 #include <float.h> /* for DBL_MAX */
29 #include <Rmath.h>
30 #include <Graphics.h>
31 #include <Print.h>
32 #include <R_ext/Boolean.h>
33
34 /* filled contours and perspective plots were originally here,
35 now in ../library/graphics/src/plot3d.c .
36 */
37
38 #include "contour-common.h"
39
40 #define CONTOUR_LIST_STEP 100
41 #define CONTOUR_LIST_LEVEL 0
42 #define CONTOUR_LIST_X 1
43 #define CONTOUR_LIST_Y 2
44
growList(SEXP oldlist)45 static SEXP growList(SEXP oldlist) {
46 int i, len;
47 SEXP templist;
48 len = LENGTH(oldlist);
49 templist = PROTECT(allocVector(VECSXP, len + CONTOUR_LIST_STEP));
50 for (i=0; i<len; i++)
51 SET_VECTOR_ELT(templist, i, VECTOR_ELT(oldlist, i));
52 UNPROTECT(1);
53 return templist;
54 }
55
56 /*
57 * Store the list of segments for a single level in the SEXP
58 * list that will be returned to the user
59 */
60 static
addContourLines(double * x,int nx,double * y,int ny,double * z,double zc,double atom,SEGP * segmentDB,int nlines,SEXP container)61 int addContourLines(double *x, int nx, double *y, int ny,
62 double *z, double zc, double atom,
63 SEGP* segmentDB, int nlines, SEXP container)
64 {
65 double xend, yend;
66 int i, ii, j, jj, ns, dir, nc;
67 SEGP seglist, seg, s, start, end;
68 SEXP ctr, level, xsxp, ysxp, names;
69 /// NB: The following is very much the same as in contour() in ../library/graphics/src/plot3d.c
70
71 /* Begin following contours. */
72 /* 1. Grab a segment */
73 /* 2. Follow its tail */
74 /* 3. Follow its head */
75 /* 4. Save the contour */
76 for (i = 0; i < nx - 1; i++)
77 for (j = 0; j < ny - 1; j++) {
78 while ((seglist = segmentDB[i + j * nx])) {
79 ii = i; jj = j;
80 start = end = seglist;
81 segmentDB[i + j * nx] = seglist->next;
82 xend = seglist->x1;
83 yend = seglist->y1;
84 while ((dir = ctr_segdir(xend, yend, x, y,
85 &ii, &jj, nx, ny))) {
86 segmentDB[ii + jj * nx]
87 = ctr_segupdate(xend, yend, dir, TRUE,/* = tail */
88 segmentDB[ii + jj * nx], &seg);
89 if (!seg) break;
90 end->next = seg;
91 end = seg;
92 xend = end->x1;
93 yend = end->y1;
94 }
95 end->next = NULL; /* <<< new for 1.2.3 */
96 ii = i; jj = j;
97 xend = seglist->x0;
98 yend = seglist->y0;
99 while ((dir = ctr_segdir(xend, yend, x, y,
100 &ii, &jj, nx, ny))) {
101 segmentDB[ii + jj * nx]
102 = ctr_segupdate(xend, yend, dir, FALSE,/* ie. head */
103 segmentDB[ii+jj*nx], &seg);
104 if (!seg) break;
105 seg->next = start;
106 start = seg;
107 xend = start->x0;
108 yend = start->y0;
109 }
110
111 /* ns := #{segments of polyline} -- need to allocate */
112 s = start;
113 ns = 0;
114 /* max_contour_segments: prevent inf.loop (shouldn't be needed) */
115 while (s && ns < max_contour_segments) {
116 ns++;
117 s = s->next;
118 }
119 if(ns == max_contour_segments)
120 warning(_("contour(): circular/long seglist -- set %s > %d?"),
121 "options(\"max.contour.segments\")", max_contour_segments);
122 /*
123 * "write" the contour locations into the list of contours
124 */
125 ctr = PROTECT(allocVector(VECSXP, 3));
126 level = PROTECT(allocVector(REALSXP, 1));
127 xsxp = PROTECT(allocVector(REALSXP, ns + 1));
128 ysxp = PROTECT(allocVector(REALSXP, ns + 1));
129 REAL(level)[0] = zc;
130 SET_VECTOR_ELT(ctr, CONTOUR_LIST_LEVEL, level);
131 s = start;
132 REAL(xsxp)[0] = s->x0;
133 REAL(ysxp)[0] = s->y0;
134 ns = 1;
135 while (s->next && ns < max_contour_segments) {
136 s = s->next;
137 REAL(xsxp)[ns] = s->x0;
138 REAL(ysxp)[ns++] = s->y0;
139 }
140 REAL(xsxp)[ns] = s->x1;
141 REAL(ysxp)[ns] = s->y1;
142 SET_VECTOR_ELT(ctr, CONTOUR_LIST_X, xsxp);
143 SET_VECTOR_ELT(ctr, CONTOUR_LIST_Y, ysxp);
144 /*
145 * Set the names attribute for the contour
146 * So that users can extract components using
147 * meaningful names
148 */
149 PROTECT(names = allocVector(STRSXP, 3));
150 SET_STRING_ELT(names, 0, mkChar("level"));
151 SET_STRING_ELT(names, 1, mkChar("x"));
152 SET_STRING_ELT(names, 2, mkChar("y"));
153 setAttrib(ctr, R_NamesSymbol, names);
154 /*
155 * We're about to add another line to the list ...
156 */
157 nlines += 1;
158 nc = LENGTH(VECTOR_ELT(container, 0));
159 if (nlines == nc)
160 /* Where does this get UNPROTECTed? */
161 SET_VECTOR_ELT(container, 0,
162 growList(VECTOR_ELT(container, 0)));
163 SET_VECTOR_ELT(VECTOR_ELT(container, 0), nlines - 1, ctr);
164 UNPROTECT(5);
165 }
166 }
167 return nlines;
168 }
169
170 /*
171 * Given nx x values, ny y values, nx*ny z values,
172 * and nl cut-values in z ...
173 * ... produce a list of contour lines:
174 * list of sub-lists
175 * sub-list = x vector, y vector, and cut-value.
176 */
GEcontourLines(double * x,int nx,double * y,int ny,double * z,double * levels,int nl)177 SEXP GEcontourLines(double *x, int nx, double *y, int ny,
178 double *z, double *levels, int nl)
179 {
180 const void *vmax;
181 int i, nlines, len;
182 double atom, zmin, zmax;
183 SEGP* segmentDB;
184 SEXP container, mainlist, templist;
185 /*
186 * "tie-breaker" values
187 */
188 zmin = DBL_MAX;
189 zmax = DBL_MIN;
190 for (i = 0; i < nx * ny; i++)
191 if (R_FINITE(z[i])) {
192 if (zmax < z[i]) zmax = z[i];
193 if (zmin > z[i]) zmin = z[i];
194 }
195
196 if (zmin >= zmax) {
197 if (zmin == zmax)
198 warning(_("all z values are equal"));
199 else
200 warning(_("all z values are NA"));
201 return R_NilValue;
202 }
203 /* change to 1e-3, reconsidered because of PR#897
204 * but 1e-7, and even 2*DBL_EPSILON do not prevent inf.loop in contour().
205 * maybe something like 16 * DBL_EPSILON * (..).
206 * see also max_contour_segments above */
207 atom = 1e-3 * (zmax - zmin);
208 /*
209 * Create a "container" which is a list with only 1 element.
210 * The element is the list of lines that will be built up.
211 * I create the container because this allows me to PROTECT
212 * the container once here and then UNPROTECT it at the end of
213 * this function and, as long as I always work with
214 * VECTOR_ELT(container, 0) and SET_VECTOR_ELT(container, 0)
215 * in functions called from here, I don't need to worry about
216 * protectin the list that I am building up.
217 * Why bother? Because the list I am building can potentially
218 * grow and it's awkward to get the PROTECTs/UNPROTECTs right
219 * when you're in a loop and growing a list.
220 */
221 container = PROTECT(allocVector(VECSXP, 1));
222 /*
223 * Create "large" list (will trim excess at the end if necesary)
224 */
225 SET_VECTOR_ELT(container, 0, allocVector(VECSXP, CONTOUR_LIST_STEP));
226 nlines = 0;
227 /*
228 * Add lines for each contour level
229 */
230 for (i = 0; i < nl; i++) {
231 /*
232 * The vmaxget/set is to manage the memory that gets
233 * R_alloc'ed in the creation of the segmentDB structure
234 */
235 vmax = vmaxget();
236 /*
237 * Generate a segment database
238 */
239 segmentDB = contourLines(x, nx, y, ny, z, levels[i], atom);
240 /*
241 * Add lines to the list based on the segment database
242 */
243 nlines = addContourLines(x, nx, y, ny, z, levels[i],
244 atom, segmentDB, nlines,
245 container);
246 vmaxset(vmax);
247 }
248 /*
249 * Trim the list of lines to the appropriate length.
250 */
251 len = LENGTH(VECTOR_ELT(container, 0));
252 if (nlines < len) {
253 mainlist = VECTOR_ELT(container, 0);
254 templist = PROTECT(allocVector(VECSXP, nlines));
255 for (i=0; i<nlines; i++)
256 SET_VECTOR_ELT(templist, i, VECTOR_ELT(mainlist, i));
257 mainlist = templist;
258 UNPROTECT(1); /* UNPROTECT templist */
259 } else
260 mainlist = VECTOR_ELT(container, 0);
261 UNPROTECT(1); /* UNPROTECT container */
262 return mainlist;
263 }
264
265 /* This is for contourLines() in package grDevices */
do_contourLines(SEXP call,SEXP op,SEXP args,SEXP env)266 SEXP do_contourLines(SEXP call, SEXP op, SEXP args, SEXP env)
267 {
268 SEXP c, x, y, z;
269 int nx, ny, nc;
270
271 x = PROTECT(coerceVector(CAR(args), REALSXP));
272 nx = LENGTH(x);
273 args = CDR(args);
274
275 y = PROTECT(coerceVector(CAR(args), REALSXP));
276 ny = LENGTH(y);
277 args = CDR(args);
278
279 z = PROTECT(coerceVector(CAR(args), REALSXP));
280 args = CDR(args);
281
282 /* levels */
283 c = PROTECT(coerceVector(CAR(args), REALSXP));
284 nc = LENGTH(c);
285 args = CDR(args);
286
287 SEXP res = GEcontourLines(REAL(x), nx, REAL(y), ny, REAL(z), REAL(c), nc);
288 UNPROTECT(4);
289 return res;
290 }
291