src/main/contour-common.h - R - Git at Google

 /*
  *  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;
 }
	/*
	* 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(nxny, 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;
	}