The R Project SVN R

/*
 *  R : A Computer Langage for Statistical Data Analysis
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 *
 *  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, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include "Defn.h"
#include "Mathlib.h"
#include "Graphics.h"

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

static SEGP *SegDB;

static int intersect(double z0, double z1, double zc, double *f)
{
    if((z0 - zc) * (z1 - zc) < 0.0) {
        *f = (zc - z0) / (z1 -  z0);
        return 1;
    }
    return 0;
}

static SEGP 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 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;
}

    /* Determine the entry direction to the next cell */
    /* and update the cell indices */

static int SegDir(double xend, double yend, double *x, double *y, int *i, int *j, int nx, int ny)
{
    if(yend == y[*j]) {
        if(*j == 0) return 0;
        *j = *j - 1;
        return 3;
    }
    if(xend == x[*i]) {
        if(*i == 0) return 0;
        *i = *i - 1;
        return 4;
    }
    if(yend == y[*j+1]) {
        if(*j >= ny - 1) return 0;
        *j = *j + 1;
        return 1;
    }
    if(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 SegUpdate(double xend, double yend, int dir, int tail, SEGP seglist, SEGP* seg)
{
    if(seglist == NULL) {
        *seg = NULL;
        return NULL;
    }
    switch(dir) {
    case 1:
    case 3:
        if(yend == seglist->y0) {
            if(!tail) SwapSeg(seglist);
            *seg = seglist;
            return seglist->next;
        }
        if(yend == seglist->y1) {
            if(tail) SwapSeg(seglist);
            *seg = seglist;
            return seglist->next;
        }
        break;
    case 2:
    case 4:
        if(xend == seglist->x0) {
            if(!tail) SwapSeg(seglist);
            *seg = seglist;
            return seglist->next;
        }
        if(xend == seglist->x1) {
            if(tail) SwapSeg(seglist);
            *seg = seglist;
            return seglist->next;
        }
        break;
    }
    seglist->next = SegUpdate(xend, yend, dir, tail, seglist->next, seg);
    return seglist;
}

static void contour(SEXP x, int nx, SEXP y, int ny, SEXP z, double zc, double atom)
{
    double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];
    double xend, yend;
    int i, ii, j, jj, k, l, m, nacode, ns, ns2, dir;
    SEGP seglist, seg, s, start, end;

    for(i=0 ; i<nx-1 ; i++) {
        xl = REAL(x)[i];
        xh = REAL(x)[i+1];
        for(j=0 ; j<ny-1 ; j++) {
            yl = REAL(y)[j];
            yh = REAL(y)[j+1];
            k = i+j*nx;
            zll = REAL(z)[k];
            zhl = REAL(z)[k+1];
            zlh = REAL(z)[k+nx];
            zhh = REAL(z)[k+nx+1];
            k = 0;

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

            if(zll == zc) zll = zll + atom;
            if(zhl == zc) zhl = zhl + atom;
            if(zlh == zc) zlh = zlh + atom;
            if(zhh == zc) zhh = zhh + atom;

                /* Check for intersections with sides */

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

            switch(nacode) {
            case 15:
                if(intersect(zll, zhl, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yl; k++;
                }
                if(intersect(zll, zlh, zc, &f)) {
                    yy[k] = yl + f * (yh - yl);
                    xx[k] = xl; k++;
                }
                if(intersect(zhl, zhh, zc, &f)) {
                    yy[k] = yl + f * (yh - yl);
                    xx[k] = xh; k++;
                }
                if(intersect(zlh, zhh, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yh; k++;
                }
                break;
            case 14:
                if(intersect(zhl, zhh, zc, &f)) {
                    yy[k] = yl + f * (yh - yl);
                    xx[k] = xh; k++;
                }
                if(intersect(zlh, zhh, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yh; k++;
                }
                if(intersect(zlh, zhl, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yh + f * (yl - yh);
                    k++;
                }
                break;
            case 13:
                if(intersect(zll, zlh, zc, &f)) {
                    yy[k] = yl + f * (yh - yl);
                    xx[k] = xl; k++;
                }
                if(intersect(zlh, zhh, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yh; k++;
                }
                if(intersect(zll, zhh, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yl + f * (yh - yl);
                    k++;
                }
                break;
            case 11:
                if(intersect(zhl, zhh, zc, &f)) {
                    yy[k] = yl + f * (yh - yl);
                    xx[k] = xh; k++;
                }
                if(intersect(zll, zhl, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yl; k++;
                }
                if(intersect(zll, zhh, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yl + f * (yh - yl);
                    k++;
                }
                break;
            case 7:
                if(intersect(zll, zlh, zc, &f)) {
                    yy[k] = yl + f * (yh - yl);
                    xx[k] = xl; k++;
                }
                if(intersect(zll, zhl, zc, &f)) {
                    xx[k] = xl + f * (xh - xl);
                    yy[k] = yl; k++;
                }
                if(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 = 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 = NewSeg(xx[0], yy[0], xx[1], yy[1], seglist);
                    seglist = NewSeg(xx[2], yy[2], xx[3], yy[3], seglist);
                }
            }
            SegDB[i+j*nx] = seglist;
        }
    }

        /* The segment database is now assembled. */
        /* Begin following contours. */
        /* 1. Grab a segment */
        /* 2. Follow its tail */
        /* 3. Follow its head */
        /* 4. Draw the contour */

    for(i=0 ; i<nx-1 ; i++)
        for(j=0 ; j<ny-1 ; j++) {
            while(seglist = SegDB[i+j*nx]) {
                ii = i; jj = j;
                start = end = seglist;
                SegDB[i+j*nx] = seglist->next;
                xend = seglist->x1;
                yend = seglist->y1;
                while(dir=SegDir(xend, yend, REAL(x), REAL(y), &ii, &jj, nx, ny)) {
                    SegDB[ii+jj*nx] = SegUpdate(xend, yend, dir, 1, SegDB[ii+jj*nx], &seg);
                    if(!seg) break;
                    end->next = seg;
                    end = seg;
                    xend = end->x1;
                    yend = end->y1;
                }
                ii = i; jj = j;
                xend = seglist->x0;
                yend = seglist->y0;
                while(dir=SegDir(xend, yend, REAL(x), REAL(y), &ii, &jj, nx, ny)) {
                    SegDB[ii+jj*nx] = SegUpdate(xend, yend, dir, 0, SegDB[ii+jj*nx], &seg);
                    if(!seg) break;
                    seg->next = start;
                    start = seg;
                    xend = start->x0;
                    yend = start->y0;
                }
                s = start;
                ns = 0;
                while(s) {
                    ns++;
                    s = s->next;
                }
                if(ns > 3) ns2 = ns/2;
                else ns2 = -1;

                s = start;
                ns = 0;
                GMode(1);
                GStartPath();
                GMoveTo(XMAP(s->x0), YMAP(s->y0));
                while(s) {
                    GLineTo(XMAP(s->x1), YMAP(s->y1));
                    s = s->next;
                }
                GEndPath();
                GMode(0);
            }
        }

}


SEXP do_contour(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP c, x, y, z, col, lty;
    int i, j, nx, ny, nc, ncol, nlty;
    int ltysave, colsave;
    double atom, zmin, zmax;
    char *vmax;

    GCheckState();

    if(length(args) < 4) errorcall(call, "too few arguments\n");

    x = CAR(args);
    internalTypeCheck(call, x, REALSXP);
    nx = LENGTH(x);
    args = CDR(args);

    y = CAR(args);
    internalTypeCheck(call, y, REALSXP);
    ny = LENGTH(y);
    args = CDR(args);

    z = CAR(args);
    internalTypeCheck(call, z, REALSXP);
    args = CDR(args);

    c = CAR(args);
    internalTypeCheck(call, c, REALSXP);
    nc = LENGTH(c);
    args = CDR(args);

    PROTECT(col = FixupCol(GetPar("col", args)));
    ncol = length(col);

    PROTECT(lty = FixupLty(GetPar("lty", args)));
    nlty = length(lty);

        /* col, lwd and lty vectors here */

    if(nx < 2 || ny < 2)
        errorcall(call, "insufficient x or y values\n");

    if(nrows(z) != nx || ncols(z) != ny)
        errorcall(call, "dimension mismatch\n");


    if(nc < 1)
        errorcall(call, "no contour values\n");

    for(i=0 ; i<nx ; i++) {
        if(!FINITE(REAL(x)[i]))
            errorcall(call, "missing x values\n");
        if(i > 0 && REAL(x)[i] < REAL(x)[i-1])
            errorcall(call, "increasing x values expected\n");
    }

    for(i=0 ; i<ny ; i++) {
        if(!FINITE(REAL(y)[i]))
            errorcall(call, "missing y values\n");
        if(i > 0 && REAL(y)[i] < REAL(y)[i-1])
            errorcall(call, "increasing y values expected\n");
    }

    for(i=0 ; i<nc ; i++)
        if(!FINITE(REAL(c)[i]))
            errorcall(call, "illegal NA contour values\n");

    zmin = DBL_MAX;
    zmax = DBL_MIN;
    for(i=0 ; i<nx*ny ; i++)
        if(FINITE(REAL(z)[i])) {
            if(zmax < REAL(z)[i]) zmax =  REAL(z)[i];
            if(zmin > REAL(z)[i]) zmin =  REAL(z)[i];
        }

    if(zmin >= zmax) {
        if(zmin == zmax)
            warning("all z values are equal\n");
        else
            warning("all z values are NA\n");
        return R_NilValue;
    }

    atom = DBL_EPSILON * (zmax - zmin);

        /* Initialize the segment data base */

    SegDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));

    for(i=0 ; i<nx ; i++)
        for(j=0 ; j<ny ; j++)
            SegDB[i+j*nx] = NULL;

        /* Draw the contours -- note the heap release */

    ltysave = GP->lty;
    colsave = GP->col;
    for(i=0 ; i<nc ; i++) {
        vmax = vmaxget();
        GP->lty = INTEGER(lty)[i%nlty];
        if(GP->lty == NA_INTEGER) GP->lty = ltysave;
        GP->col = INTEGER(col)[i%ncol];
        if(GP->col == NA_INTEGER) GP->col = colsave;
        contour(x, nx, y, ny, z, REAL(c)[i], atom);
        vmaxset(vmax);
    }
    GP->lty = ltysave;
    GP->col = colsave;
    UNPROTECT(2);
    return R_NilValue;
}


SEXP do_image(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP sx, sy, sz, szlim, sc;
    double *x, *y, *z;
    unsigned *c;
    double xlow, xhigh, ylow, yhigh, zmin, zmax;
    int i, j, nx, ny, nz, ic, nc, colsave, xpdsave;

    GCheckState();

    checkArity(op,args);

    sx = CAR(args);
    internalTypeCheck(call, sx, REALSXP);
    nx = LENGTH(sx);
    args = CDR(args);

    sy = CAR(args);
    internalTypeCheck(call, sy, REALSXP);
    ny = LENGTH(sy);
    args = CDR(args);

    sz = CAR(args);
    internalTypeCheck(call, sz, REALSXP);
    nz = length(sz);
    args = CDR(args);

    szlim = CAR(args);
    internalTypeCheck(call, szlim, REALSXP);
    if(length(szlim) != 2 ||
        !FINITE(REAL(szlim)[0]) ||
        !FINITE(REAL(szlim)[1]) ||
        REAL(szlim)[0] >= REAL(szlim)[1])
            errorcall(call, "invalid z limits\n");
    zmin = REAL(szlim)[0];
    zmax = REAL(szlim)[1];
    args = CDR(args);

    PROTECT(sc = FixupCol(CAR(args)));
    nc = length(sc);

        /* Shorthand Pointers */

    x = REAL(sx);
    y = REAL(sy);
    z = REAL(sz);
    c = (unsigned*)INTEGER(sc);

        /* Check of grid coordinates */
        /* We want them to all be finite and */
        /* in strictly ascending order */

    if(nx < 2 || ny < 2) goto badxy;
    if(!FINITE(x[0])) goto badxy;
    if(!FINITE(y[0])) goto badxy;
    for(i=1 ; i<nx ; i++)
        if(!FINITE(x[i]) || x[i] <= x[i-1]) goto badxy;
    for(j=1 ; j<ny ; j++)
        if(!FINITE(y[j]) || y[j] <= y[j-1]) goto badxy;

    colsave = GP->col;
    xpdsave = GP->xpd;
    GP->xpd = 0;

    GMode(1);

    for(i=0 ; i<nx ; i++) {
        if(i == 0)
            xlow = XMAP(x[0]);
        else
            xlow = XMAP(0.5 * (x[i] + x[i-1]));
        if(i == nx-1)
            xhigh = XMAP(x[nx-1]);
        else
            xhigh = XMAP(0.5 * (x[i] + x[i+1]));

        for(j=0 ; j<ny ; j++) {
            if(FINITE(z[i+j*nx])) {
                ic = floor((nc - 1) * (z[i+j*nx]-zmin)/(zmax - zmin) + 0.5);
                if(ic >= 0 && ic < nc) {
                    if(j == 0)
                        ylow = YMAP(y[0]);
                    else
                        ylow = YMAP(0.5 * (y[j] + y[j-1]));
                    if(j == ny-1)
                        yhigh = YMAP(y[ny-1]);
                    else
                        yhigh = YMAP(0.5 * (y[j] + y[j+1]));
                    GRect(xlow, ylow, xhigh, yhigh, c[ic], NA_INTEGER);
                }
            }
        }
    }
    GMode(0);
    GP->col = colsave;
    GP->xpd = xpdsave;
    R_Visible = 0;
    UNPROTECT(1);
    return R_NilValue;

badxy:
    errorcall(call, "invalid x / y limits\n");
}