| /* |
| * R : A Computer Language for Statistical Data Analysis |
| * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka |
| * Copyright (C) 1997--2014 The R Core Team |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| */ |
| |
| |
| #ifdef HAVE_CONFIG_H |
| #include <config.h> |
| #endif |
| |
| #include <Defn.h> |
| #include <Internal.h> |
| #include <float.h> /* for DBL_MAX */ |
| #include <Rmath.h> |
| #include <Graphics.h> |
| #include <Print.h> |
| #include <R_ext/Boolean.h> |
| |
| /* filled contours and perspective plots were originally here, |
| now in graphics/src/plot3d.c . |
| */ |
| |
| #include "contour-common.h" |
| |
| #define CONTOUR_LIST_STEP 100 |
| #define CONTOUR_LIST_LEVEL 0 |
| #define CONTOUR_LIST_X 1 |
| #define CONTOUR_LIST_Y 2 |
| |
| static SEXP growList(SEXP oldlist) { |
| int i, len; |
| SEXP templist; |
| len = LENGTH(oldlist); |
| templist = PROTECT(allocVector(VECSXP, len + CONTOUR_LIST_STEP)); |
| for (i=0; i<len; i++) |
| SET_VECTOR_ELT(templist, i, VECTOR_ELT(oldlist, i)); |
| UNPROTECT(1); |
| return templist; |
| } |
| |
| /* |
| * Store the list of segments for a single level in the SEXP |
| * list that will be returned to the user |
| */ |
| static |
| int addContourLines(double *x, int nx, double *y, int ny, |
| double *z, double zc, double atom, |
| SEGP* segmentDB, int nlines, SEXP container) |
| { |
| double xend, yend; |
| int i, ii, j, jj, ns, dir, nc; |
| SEGP seglist, seg, s, start, end; |
| SEXP ctr, level, xsxp, ysxp, names; |
| /* Begin following contours. */ |
| /* 1. Grab a segment */ |
| /* 2. Follow its tail */ |
| /* 3. Follow its head */ |
| /* 4. Save the contour */ |
| for (i = 0; i < nx - 1; i++) |
| for (j = 0; j < ny - 1; j++) { |
| while ((seglist = segmentDB[i + j * nx])) { |
| ii = i; jj = j; |
| start = end = seglist; |
| segmentDB[i + j * nx] = seglist->next; |
| xend = seglist->x1; |
| yend = seglist->y1; |
| while ((dir = ctr_segdir(xend, yend, x, y, |
| &ii, &jj, nx, ny))) { |
| segmentDB[ii + jj * nx] |
| = ctr_segupdate(xend, yend, dir, TRUE,/* = tail */ |
| segmentDB[ii + jj * nx], &seg); |
| if (!seg) break; |
| end->next = seg; |
| end = seg; |
| xend = end->x1; |
| yend = end->y1; |
| } |
| end->next = NULL; /* <<< new for 1.2.3 */ |
| ii = i; jj = j; |
| xend = seglist->x0; |
| yend = seglist->y0; |
| while ((dir = ctr_segdir(xend, yend, x, y, |
| &ii, &jj, nx, ny))) { |
| segmentDB[ii + jj * nx] |
| = ctr_segupdate(xend, yend, dir, FALSE,/* ie. head */ |
| segmentDB[ii+jj*nx], &seg); |
| if (!seg) break; |
| seg->next = start; |
| start = seg; |
| xend = start->x0; |
| yend = start->y0; |
| } |
| |
| /* ns := #{segments of polyline} -- need to allocate */ |
| s = start; |
| ns = 0; |
| /* max_contour_segments: prevent inf.loop (shouldn't be needed) */ |
| while (s && ns < max_contour_segments) { |
| ns++; |
| s = s->next; |
| } |
| if(ns == max_contour_segments) |
| warning(_("contour(): circular/long seglist -- set %s > %d?"), |
| "options(\"max.contour.segments\")", max_contour_segments); |
| /* |
| * "write" the contour locations into the list of contours |
| */ |
| ctr = PROTECT(allocVector(VECSXP, 3)); |
| level = PROTECT(allocVector(REALSXP, 1)); |
| xsxp = PROTECT(allocVector(REALSXP, ns + 1)); |
| ysxp = PROTECT(allocVector(REALSXP, ns + 1)); |
| REAL(level)[0] = zc; |
| SET_VECTOR_ELT(ctr, CONTOUR_LIST_LEVEL, level); |
| s = start; |
| REAL(xsxp)[0] = s->x0; |
| REAL(ysxp)[0] = s->y0; |
| ns = 1; |
| while (s->next && ns < max_contour_segments) { |
| s = s->next; |
| REAL(xsxp)[ns] = s->x0; |
| REAL(ysxp)[ns++] = s->y0; |
| } |
| REAL(xsxp)[ns] = s->x1; |
| REAL(ysxp)[ns] = s->y1; |
| SET_VECTOR_ELT(ctr, CONTOUR_LIST_X, xsxp); |
| SET_VECTOR_ELT(ctr, CONTOUR_LIST_Y, ysxp); |
| /* |
| * Set the names attribute for the contour |
| * So that users can extract components using |
| * meaningful names |
| */ |
| PROTECT(names = allocVector(STRSXP, 3)); |
| SET_STRING_ELT(names, 0, mkChar("level")); |
| SET_STRING_ELT(names, 1, mkChar("x")); |
| SET_STRING_ELT(names, 2, mkChar("y")); |
| setAttrib(ctr, R_NamesSymbol, names); |
| /* |
| * We're about to add another line to the list ... |
| */ |
| nlines += 1; |
| nc = LENGTH(VECTOR_ELT(container, 0)); |
| if (nlines == nc) |
| /* Where does this get UNPROTECTed? */ |
| SET_VECTOR_ELT(container, 0, |
| growList(VECTOR_ELT(container, 0))); |
| SET_VECTOR_ELT(VECTOR_ELT(container, 0), nlines - 1, ctr); |
| UNPROTECT(5); |
| } |
| } |
| return nlines; |
| } |
| |
| /* |
| * Given nx x values, ny y values, nx*ny z values, |
| * and nl cut-values in z ... |
| * ... produce a list of contour lines: |
| * list of sub-lists |
| * sub-list = x vector, y vector, and cut-value. |
| */ |
| SEXP GEcontourLines(double *x, int nx, double *y, int ny, |
| double *z, double *levels, int nl) |
| { |
| const void *vmax; |
| int i, nlines, len; |
| double atom, zmin, zmax; |
| SEGP* segmentDB; |
| SEXP container, mainlist, templist; |
| /* |
| * "tie-breaker" values |
| */ |
| zmin = DBL_MAX; |
| zmax = DBL_MIN; |
| for (i = 0; i < nx * ny; i++) |
| if (R_FINITE(z[i])) { |
| if (zmax < z[i]) zmax = z[i]; |
| if (zmin > z[i]) zmin = z[i]; |
| } |
| |
| if (zmin >= zmax) { |
| if (zmin == zmax) |
| warning(_("all z values are equal")); |
| else |
| warning(_("all z values are NA")); |
| return R_NilValue; |
| } |
| /* change to 1e-3, reconsidered because of PR#897 |
| * but 1e-7, and even 2*DBL_EPSILON do not prevent inf.loop in contour(). |
| * maybe something like 16 * DBL_EPSILON * (..). |
| * see also max_contour_segments above */ |
| atom = 1e-3 * (zmax - zmin); |
| /* |
| * Create a "container" which is a list with only 1 element. |
| * The element is the list of lines that will be built up. |
| * I create the container because this allows me to PROTECT |
| * the container once here and then UNPROTECT it at the end of |
| * this function and, as long as I always work with |
| * VECTOR_ELT(container, 0) and SET_VECTOR_ELT(container, 0) |
| * in functions called from here, I don't need to worry about |
| * protectin the list that I am building up. |
| * Why bother? Because the list I am building can potentially |
| * grow and it's awkward to get the PROTECTs/UNPROTECTs right |
| * when you're in a loop and growing a list. |
| */ |
| container = PROTECT(allocVector(VECSXP, 1)); |
| /* |
| * Create "large" list (will trim excess at the end if necesary) |
| */ |
| SET_VECTOR_ELT(container, 0, allocVector(VECSXP, CONTOUR_LIST_STEP)); |
| nlines = 0; |
| /* |
| * Add lines for each contour level |
| */ |
| for (i = 0; i < nl; i++) { |
| /* |
| * The vmaxget/set is to manage the memory that gets |
| * R_alloc'ed in the creation of the segmentDB structure |
| */ |
| vmax = vmaxget(); |
| /* |
| * Generate a segment database |
| */ |
| segmentDB = contourLines(x, nx, y, ny, z, levels[i], atom); |
| /* |
| * Add lines to the list based on the segment database |
| */ |
| nlines = addContourLines(x, nx, y, ny, z, levels[i], |
| atom, segmentDB, nlines, |
| container); |
| vmaxset(vmax); |
| } |
| /* |
| * Trim the list of lines to the appropriate length. |
| */ |
| len = LENGTH(VECTOR_ELT(container, 0)); |
| if (nlines < len) { |
| mainlist = VECTOR_ELT(container, 0); |
| templist = PROTECT(allocVector(VECSXP, nlines)); |
| for (i=0; i<nlines; i++) |
| SET_VECTOR_ELT(templist, i, VECTOR_ELT(mainlist, i)); |
| mainlist = templist; |
| UNPROTECT(1); /* UNPROTECT templist */ |
| } else |
| mainlist = VECTOR_ELT(container, 0); |
| UNPROTECT(1); /* UNPROTECT container */ |
| return mainlist; |
| } |
| |
| /* This is for contourLines() in package grDevices */ |
| SEXP do_contourLines(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP c, x, y, z; |
| int nx, ny, nc; |
| |
| x = PROTECT(coerceVector(CAR(args), REALSXP)); |
| nx = LENGTH(x); |
| args = CDR(args); |
| |
| y = PROTECT(coerceVector(CAR(args), REALSXP)); |
| ny = LENGTH(y); |
| args = CDR(args); |
| |
| z = PROTECT(coerceVector(CAR(args), REALSXP)); |
| args = CDR(args); |
| |
| /* levels */ |
| c = PROTECT(coerceVector(CAR(args), REALSXP)); |
| nc = LENGTH(c); |
| args = CDR(args); |
| |
| SEXP res = GEcontourLines(REAL(x), nx, REAL(y), ny, REAL(z), REAL(c), nc); |
| UNPROTECT(4); |
| return res; |
| } |