The R Project SVN R

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 *  Copyright (C) 1997--2013  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/
 */


/* Stuff for labels on contour plots
   Originally written by Nicholas Hildreth
   Adapted by Paul Murrell
*/

/* Included by src/main/plot3d.c and src/library/graphics/src/plot3d  */

    /*  C o n t o u r   P l o t t i n g  */

typedef struct SEG {
    struct SEG *next;
    double x0;
    double y0;
    double x1;
    double y1;
} SEG, *SEGP;


static int ctr_intersect(double z0, double z1, double zc, double *f)
{
/*  Old test was  ((z0 - zc) * (z1 - zc) < 0.0), but rounding led to inconsistencies
    in PR#15454 */
    if ( (z0 < zc) != (z1 < zc) && z0 != zc && z1 != zc ) {
    *f = (zc - z0) / (z1 -  z0);
    return 1;
    }
    return 0;
}

static SEGP ctr_newseg(double x0, double y0, double x1, double y1, SEGP prev)
{
    SEGP seg = (SEGP)R_alloc(1, sizeof(SEG));
    seg->x0 = x0;
    seg->y0 = y0;
    seg->x1 = x1;
    seg->y1 = y1;
    seg->next = prev;
    return seg;
}

static void ctr_swapseg(SEGP seg)
{
    double x, y;
    x = seg->x0;
    y = seg->y0;
    seg->x0 = seg->x1;
    seg->y0 = seg->y1;
    seg->x1 = x;
    seg->y1 = y;
}

    /* ctr_segdir(): Determine the entry direction to the next cell */
    /* and update the cell indices */

#define XMATCH(x0,x1) (fabs(x0-x1) == 0)
#define YMATCH(y0,y1) (fabs(y0-y1) == 0)

static int ctr_segdir(double xend, double yend, double *x, double *y,
              int *i, int *j, int nx, int ny)
{
    if (YMATCH(yend, y[*j])) {
    if (*j == 0)
        return 0;
    *j = *j - 1;
    return 3;
    }
    if (XMATCH(xend, x[*i])) {
    if (*i == 0)
        return 0;
    *i = *i - 1;
    return 4;
    }
    if (YMATCH(yend, y[*j + 1])) {
    if (*j >= ny - 1)
        return 0;
    *j = *j + 1;
    return 1;
    }
    if (XMATCH(xend, x[*i + 1])) {
    if (*i >= nx - 1)
        return 0;
    *i = *i + 1;
    return 2;
    }
    return 0;
}

/* Search seglist for a segment with endpoint (xend, yend). */
/* The cell entry direction is dir, and if tail=1/0 we are */
/* building the tail/head of a contour.  The matching segment */
/* is pointed to by seg and the updated segment list (with */
/* the matched segment stripped) is returned by the funtion. */

static SEGP ctr_segupdate(double xend, double yend, int dir, Rboolean tail,
              SEGP seglist, SEGP* seg)
{
    if (seglist == NULL) {
    *seg = NULL;
    return NULL;
    }
    switch (dir) {
    case 1:
    case 3:
    if (YMATCH(yend,seglist->y0)) {
        if (!tail)
        ctr_swapseg(seglist);
        *seg = seglist;
        return seglist->next;
    }
    if (YMATCH(yend,seglist->y1)) {
        if (tail)
        ctr_swapseg(seglist);
        *seg = seglist;
        return seglist->next;
    }
    break;
    case 2:
    case 4:
    if (XMATCH(xend,seglist->x0)) {
        if (!tail)
        ctr_swapseg(seglist);
        *seg = seglist;
        return seglist->next;
    }
    if (XMATCH(xend,seglist->x1)) {
        if (tail)
        ctr_swapseg(seglist);
        *seg = seglist;
        return seglist->next;
    }
    break;
    }
    seglist->next = ctr_segupdate(xend, yend, dir, tail, seglist->next, seg);
    return seglist;
}



/*
 * Generate a list of segments for a single level
 *
 * NB this R_allocs its return value, so callers need to manage R_alloc stack.
 */
static SEGP* contourLines(double *x, int nx, double *y, int ny,
             double *z, double zc, double atom)
{
    double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];
    int i, j, k, l, m, nacode;
    SEGP seglist;
    SEGP *segmentDB;
    /* Initialize the segment data base */
    /* Note we must be careful about resetting */
    /* the top of the stack, otherwise we run out of */
    /* memory after a sequence of displaylist replays */
    /*
     * This reset is done out in GEcontourLines
     */
    segmentDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));
    for (i = 0; i < nx; i++)
    for (j = 0; j < ny; j++)
        segmentDB[i + j * nx] = NULL;
    for (i = 0; i < nx - 1; i++) {
    xl = x[i];
    xh = x[i + 1];
    for (j = 0; j < ny - 1; j++) {
        yl = y[j];
        yh = y[j + 1];
        k = i + j * nx;
        zll = z[k];
        zhl = z[k + 1];
        zlh = z[k + nx];
        zhh = z[k + nx + 1];

        /* If the value at a corner is exactly equal to a contour level,
         * change that value by a tiny amount */

        if (zll == zc) zll += atom;
        if (zhl == zc) zhl += atom;
        if (zlh == zc) zlh += atom;
        if (zhh == zc) zhh += atom;
#ifdef DEBUG_contour
        /* Haven't seen this happening (MM): */
        if (zll == zc) REprintf(" [%d,%d] ll: %g\n",i,j, zll);
        if (zhl == zc) REprintf(" [%d,%d] hl: %g\n",i,j, zhl);
        if (zlh == zc) REprintf(" [%d,%d] lh: %g\n",i,j, zlh);
        if (zhh == zc) REprintf(" [%d,%d] hh: %g\n",i,j, zhh);
#endif
        /* Check for intersections with sides */

        nacode = 0;
        if (R_FINITE(zll)) nacode += 1;
        if (R_FINITE(zhl)) nacode += 2;
        if (R_FINITE(zlh)) nacode += 4;
        if (R_FINITE(zhh)) nacode += 8;

        k = 0;
        switch (nacode) {
        case 15:
        if (ctr_intersect(zll, zhl, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yl; k++;
        }
        if (ctr_intersect(zll, zlh, zc, &f)) {
            yy[k] = yl + f * (yh - yl);
            xx[k] = xl; k++;
        }
        if (ctr_intersect(zhl, zhh, zc, &f)) {
            yy[k] = yl + f * (yh - yl);
            xx[k] = xh; k++;
        }
        if (ctr_intersect(zlh, zhh, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yh; k++;
        }
        break;
        case 14:
        if (ctr_intersect(zhl, zhh, zc, &f)) {
            yy[k] = yl + f * (yh - yl);
            xx[k] = xh; k++;
        }
        if (ctr_intersect(zlh, zhh, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yh; k++;
        }
        if (ctr_intersect(zlh, zhl, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yh + f * (yl - yh);
            k++;
        }
        break;
        case 13:
        if (ctr_intersect(zll, zlh, zc, &f)) {
            yy[k] = yl + f * (yh - yl);
            xx[k] = xl; k++;
        }
        if (ctr_intersect(zlh, zhh, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yh; k++;
        }
        if (ctr_intersect(zll, zhh, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yl + f * (yh - yl);
            k++;
        }
        break;
        case 11:
        if (ctr_intersect(zhl, zhh, zc, &f)) {
            yy[k] = yl + f * (yh - yl);
            xx[k] = xh; k++;
        }
        if (ctr_intersect(zll, zhl, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yl; k++;
        }
        if (ctr_intersect(zll, zhh, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yl + f * (yh - yl);
            k++;
        }
        break;
        case 7:
        if (ctr_intersect(zll, zlh, zc, &f)) {
            yy[k] = yl + f * (yh - yl);
            xx[k] = xl; k++;
        }
        if (ctr_intersect(zll, zhl, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yl; k++;
        }
        if (ctr_intersect(zlh, zhl, zc, &f)) {
            xx[k] = xl + f * (xh - xl);
            yy[k] = yh + f * (yl - yh);
            k++;
        }
        break;
        }

        /* We now have k(=2,4) endpoints */
        /* Decide which to join */

        seglist = NULL;

        if (k > 0) {
        if (k == 2) {
            seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
        }
        else if (k == 4) {
            for (k = 3; k >= 1; k--) {
            m = k;
            xl = xx[k];
            for (l = 0; l < k; l++) {
                if (xx[l] > xl) {
                xl = xx[l];
                m = l;
                }
            }
            if (m != k) {
                xl = xx[k];
                yl = yy[k];
                xx[k] = xx[m];
                yy[k] = yy[m];
                xx[m] = xl;
                yy[m] = yl;
            }
            }
            seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
            seglist = ctr_newseg(xx[2], yy[2], xx[3], yy[3], seglist);
        }
        else error("k = %d, should be 2 or 4", k);
        }
        segmentDB[i + j * nx] = seglist;
    }
    }
    return segmentDB;
}