home *** CD-ROM | disk | FTP | other *** search
- #ifndef lint
- static char *RCSid = "$Id: util3d.c,v 1.8 1995/12/11 23:14:08 drd Exp $";
- #endif
-
-
- /* GNUPLOT - util3d.c */
- /*
- * Copyright (C) 1986 - 1993 Thomas Williams, Colin Kelley
- *
- * Permission to use, copy, and distribute this software and its
- * documentation for any purpose with or without fee is hereby granted,
- * provided that the above copyright notice appear in all copies and
- * that both that copyright notice and this permission notice appear
- * in supporting documentation.
- *
- * Permission to modify the software is granted, but not the right to
- * distribute the modified code. Modifications are to be distributed
- * as patches to released version.
- *
- * This software is provided "as is" without express or implied warranty.
- *
- *
- * AUTHORS
- *
- * Original Software:
- * Gershon Elber and many others.
- *
- * 19 September 1992 Lawrence Crowl (crowl@cs.orst.edu)
- * Added user-specified bases for log scaling.
- *
- * 3.6 - split graph3d.c into graph3d.c (graph),
- * util3d.c (intersections, etc)
- * hidden3d.c (hidden-line removal code)
- *
- * There is a mailing list for gnuplot users. Note, however, that the
- * newsgroup
- * comp.graphics.gnuplot
- * is identical to the mailing list (they
- * both carry the same set of messages). We prefer that you read the
- * messages through that newsgroup, to subscribing to the mailing list.
- * (If you can read that newsgroup, and are already on the mailing list,
- * please send a message info-gnuplot-request@dartmouth.edu, asking to be
- * removed from the mailing list.)
- *
- * The address for mailing to list members is
- * info-gnuplot@dartmouth.edu
- * and for mailing administrative requests is
- * info-gnuplot-request@dartmouth.edu
- * The mailing list for bug reports is
- * bug-gnuplot@dartmouth.edu
- * The list of those interested in beta-test versions is
- * info-gnuplot-beta@dartmouth.edu
- */
-
- #include <math.h>
- #if !defined(sequent) && !defined(apollo) && !defined(alliant)
- #include <limits.h>
- #endif
- #include "plot.h"
- #include "setshow.h"
-
- extern int xleft,xright,ybot,ytop;
- extern int hidden_active, hidden_no_update;
-
- /* ACCESS THINGS THAT OUGHT TO BE HIDDEN IN hidden3d.c - perhaps we
- * can move the relevant code into hidden3d.c sometime
- */
-
- extern long int xmin_hl,xmax_hl;
- /* These arrays are used to keep track of the minimum and maximum y values used
- for each X value. These are only used for drawing the individual boxes that
- make up the 3d figure. After each box is drawn, the information is copied
- to the bitmap. */
- extern short int *ymin_hl, *ymax_hl;
- extern int xfact, yfact;
- #define XREDUCE(X) ((X)/xfact)
- #define YREDUCE(Y) ((Y)/yfact)
- /* Bitmap of the screen. The array for each x value is malloc-ed as needed */
- extern short int **pnt;
- #define IFSET(X,Y) (pnt[X] == 0 ? 0 : (((pnt[X])[(Y)>>4] >> ((Y) & 0xf)) & 0x01))
-
- extern int suppressMove;
-
-
- extern double min_array[], max_array[];
- extern int auto_array[], log_array[];
- extern double base_array[], log_base_array[];
-
- /* for convenience while converting to use these arrays */
- #define x_min3d min_array[FIRST_X_AXIS]
- #define x_max3d max_array[FIRST_X_AXIS]
- #define y_min3d min_array[FIRST_Y_AXIS]
- #define y_max3d max_array[FIRST_Y_AXIS]
- #define z_min3d min_array[FIRST_Z_AXIS]
- #define z_max3d max_array[FIRST_Z_AXIS]
- #define min3d_z min_array[FIRST_Z_AXIS]
- #define max3d_z max_array[FIRST_Z_AXIS]
-
- /* both min() and MIN() are predefined by some compilers :-( */
-
- #define GPMIN(a,b) ( (a) < (b) ? (a) : (b) )
- #define GPMAX(a,b) ( (a) > (b) ? (a) : (b) )
-
- #define inrange(z,min,max) ((min<max) ? ((z>=min)&&(z<=max)) : ((z>=max)&&(z<=min)) )
-
-
- typedef double transform_matrix[4][4];
-
- void mat_unit(mat)
- transform_matrix mat;
- {
- int i, j;
-
- for (i = 0; i < 4; i++) for (j = 0; j < 4; j++)
- if (i == j)
- mat[i][j] = 1.0;
- else
- mat[i][j] = 0.0;
- }
-
- void mat_trans(tx, ty, tz, mat)
- double tx, ty, tz;
- transform_matrix mat;
- {
- mat_unit(mat); /* Make it unit matrix. */
- mat[3][0] = tx;
- mat[3][1] = ty;
- mat[3][2] = tz;
- }
-
- void mat_scale(sx, sy, sz, mat)
- double sx, sy, sz;
- transform_matrix mat;
- {
- mat_unit(mat); /* Make it unit matrix. */
- mat[0][0] = sx;
- mat[1][1] = sy;
- mat[2][2] = sz;
- }
-
- void mat_rot_x(teta, mat)
- double teta;
- transform_matrix mat;
- {
- double cos_teta, sin_teta;
-
- teta *= Pi / 180.0;
- cos_teta = cos(teta);
- sin_teta = sin(teta);
-
- mat_unit(mat); /* Make it unit matrix. */
- mat[1][1] = cos_teta;
- mat[1][2] = -sin_teta;
- mat[2][1] = sin_teta;
- mat[2][2] = cos_teta;
- }
-
- void mat_rot_y(teta, mat)
- double teta;
- transform_matrix mat;
- {
- double cos_teta, sin_teta;
-
- teta *= Pi / 180.0;
- cos_teta = cos(teta);
- sin_teta = sin(teta);
-
- mat_unit(mat); /* Make it unit matrix. */
- mat[0][0] = cos_teta;
- mat[0][2] = -sin_teta;
- mat[2][0] = sin_teta;
- mat[2][2] = cos_teta;
- }
-
- void mat_rot_z(teta, mat)
- double teta;
- transform_matrix mat;
- {
- double cos_teta, sin_teta;
-
- teta *= Pi / 180.0;
- cos_teta = cos(teta);
- sin_teta = sin(teta);
-
- mat_unit(mat); /* Make it unit matrix. */
- mat[0][0] = cos_teta;
- mat[0][1] = -sin_teta;
- mat[1][0] = sin_teta;
- mat[1][1] = cos_teta;
- }
-
- /* Multiply two transform_matrix. Result can be one of two operands. */
- void mat_mult(mat_res, mat1, mat2)
- transform_matrix mat_res, mat1, mat2;
- {
- int i, j, k;
- transform_matrix mat_res_temp;
-
- for (i = 0; i < 4; i++) for (j = 0; j < 4; j++) {
- mat_res_temp[i][j] = 0;
- for (k = 0; k < 4; k++) mat_res_temp[i][j] += mat1[i][k] * mat2[k][j];
- }
- for (i = 0; i < 4; i++) for (j = 0; j < 4; j++)
- mat_res[i][j] = mat_res_temp[i][j];
- }
-
-
- /* Test a single point to be within the xleft,xright,ybot,ytop bbox.
- * Sets the returned integers 4 l.s.b. as follows:
- * bit 0 if to the left of xleft.
- * bit 1 if to the right of xright.
- * bit 2 if above of ytop.
- * bit 3 if below of ybot.
- * 0 is returned if inside.
- */
- int clip_point(x, y)
- unsigned int x, y;
- {
- int ret_val = 0;
-
- if (x < xleft) ret_val |= 0x01;
- if (x > xright) ret_val |= 0x02;
- if (y < ybot) ret_val |= 0x04;
- if (y > ytop) ret_val |= 0x08;
-
- return ret_val;
- }
-
- /* Clip the given line to drawing coords defined as xleft,xright,ybot,ytop.
- * This routine uses the cohen & sutherland bit mapping for fast clipping -
- * see "Principles of Interactive Computer Graphics" Newman & Sproull page 65.
- */
- void draw_clip_line(x1, y1, x2, y2)
- unsigned int x1, y1, x2, y2;
- {
- int x, y, dx, dy, x_intr[4], y_intr[4], count, pos1, pos2;
- register struct termentry *t = term;
-
- #if defined(ATARI) || defined(MTOS)
- if(x1<0||x2<0||y1<0||y2<0) return; /* temp bug fix */
- #endif
-
- pos1 = clip_point(x1, y1);
- pos2 = clip_point(x2, y2);
- if (pos1 || pos2) {
- if (pos1 & pos2) return; /* segment is totally out. */
-
- /* Here part of the segment MAY be inside. test the intersection
- * of this segment with the 4 boundaries for hopefully 2 intersections
- * in. If none are found segment is totaly out.
- * Under rare circumstances there may be up to 4 intersections (e.g.
- * when the line passes directly through at least one corner). In
- * this case it is sufficient to take any 2 intersections (e.g. the
- * first two found).
- */
- count = 0;
- dx = x2 - x1;
- dy = y2 - y1;
-
- /* Find intersections with the x parallel bbox lines: */
- if (dy != 0) {
- x = (ybot - y2) * dx / dy + x2; /* Test for ybot boundary. */
- if (x >= xleft && x <= xright) {
- x_intr[count] = x;
- y_intr[count++] = ybot;
- }
- x = (ytop - y2) * dx / dy + x2; /* Test for ytop boundary. */
- if (x >= xleft && x <= xright) {
- x_intr[count] = x;
- y_intr[count++] = ytop;
- }
- }
-
- /* Find intersections with the y parallel bbox lines: */
- if (dx != 0) {
- y = (xleft - x2) * dy / dx + y2; /* Test for xleft boundary. */
- if (y >= ybot && y <= ytop) {
- x_intr[count] = xleft;
- y_intr[count++] = y;
- }
- y = (xright - x2) * dy / dx + y2; /* Test for xright boundary. */
- if (y >= ybot && y <= ytop) {
- x_intr[count] = xright;
- y_intr[count++] = y;
- }
- }
-
- if (count >= 2) {
- int x_max, x_min, y_max, y_min;
-
- x_min = GPMIN(x1, x2);
- x_max = GPMAX(x1, x2);
- y_min = GPMIN(y1, y2);
- y_max = GPMAX(y1, y2);
-
- if (pos1 && pos2) { /* Both were out - update both */
- x1 = x_intr[0];
- y1 = y_intr[0];
- x2 = x_intr[1];
- y2 = y_intr[1];
- }
- else if (pos1) { /* Only x1/y1 was out - update only it */
- if (dx * (x2 - x_intr[0]) + dy * (y2 - y_intr[0]) > 0) {
- x1 = x_intr[0];
- y1 = y_intr[0];
- }
- else {
- x1 = x_intr[1];
- y1 = y_intr[1];
- }
- }
- else { /* Only x2/y2 was out - update only it */
- if (dx * (x_intr[0] - x1) + dy * (y_intr[0] - y1) > 0) {
- x2 = x_intr[0];
- y2 = y_intr[0];
- }
- else {
- x2 = x_intr[1];
- y2 = y_intr[1];
- }
- }
-
- if (x1 < x_min || x1 > x_max ||
- x2 < x_min || x2 > x_max ||
- y1 < y_min || y1 > y_max ||
- y2 < y_min || y2 > y_max) return;
- }
- else
- return;
- }
-
- /* I'D REALLY LIKE TO MOVE THIS INTO hidden3d.c AS A FN - it shouldn't
- * be accessing secret variables like this
- */
-
- #ifndef LITE
- if(hidden3d && hidden_active && draw_surface)
- {
- char flag;
- register int xv, yv, errx, erry, err;
- register int xvr, yvr;
- int xve, yve;
- register int dy, nstep=0, dyr;
- int i;
- if (x1 > x2){
- xvr = x2;
- yvr = y2;
- xve = x1;
- yve = y1;
- } else {
- xvr = x1;
- yvr = y1;
- xve = x2;
- yve = y2;
- };
- errx = XREDUCE(xve) - XREDUCE(xvr);
- erry = YREDUCE(yve) - YREDUCE(yvr);
- dy = (erry > 0 ? 1 : -1);
- dyr = dy*yfact;
- switch (dy){
- case 1:
- nstep = errx + erry;
- errx = -errx;
- break;
- case -1:
- nstep = errx - erry;
- errx = -errx;
- erry = -erry;
- break;
- };
- err = errx + erry;
- errx <<= 1;
- erry <<= 1;
- xv = XREDUCE(xvr) - XREDUCE(xleft);
- yv = YREDUCE(yvr) - YREDUCE(ybot);
- (*t->move)(xvr,yvr);
- if( !IFSET(xv,yv) ) flag = 0;
- else flag = 1;
- if(!hidden_no_update){ /* Check first point */
- if (xv < xmin_hl) xmin_hl = xv;
- if (xv > xmax_hl) xmax_hl = xv;
- if (yv > ymax_hl[xv]) ymax_hl[xv] = yv;
- if (yv < ymin_hl[xv]) ymin_hl[xv] = yv;
- };
- for (i=0;i<nstep;i++){
- if (err < 0){
- xv ++;
- xvr += xfact;
- err += erry;
- } else {
- yv += dy;
- yvr += dyr;
- err += errx;
- };
- if( !IFSET(xv,yv)){
- if(flag != 0) {(*t->move)(xvr,yvr); flag = 0;};
- } else {
- if(flag == 0) {(*t->vector)(xvr,yvr); flag = 1;};
- };
- if(!hidden_no_update){
- if (xv < xmin_hl) xmin_hl = xv;
- if (xv > xmax_hl) xmax_hl = xv;
- if (yv > ymax_hl[xv]) ymax_hl[xv] = yv;
- if (yv < ymin_hl[xv]) ymin_hl[xv] = yv;
- };
- };
- if (flag == 0) (*t->vector)(xve, yve);
- return;
- };
- #endif /* not LITE */
- if(!suppressMove) (*t->move)(x1,y1);
- (*t->vector)(x2,y2);
- }
-
- /* Two routine to emulate move/vector sequence using line drawing routine. */
- static unsigned int move_pos_x, move_pos_y;
-
- void clip_move(x,y)
- unsigned int x,y;
- {
- move_pos_x = x;
- move_pos_y = y;
- }
-
- void clip_vector(x,y)
- unsigned int x,y;
- {
- draw_clip_line(move_pos_x,move_pos_y, x, y);
- move_pos_x = x;
- move_pos_y = y;
- }
-
- /* And text clipping routine. */
- void clip_put_text(x, y, str)
- unsigned int x,y;
- char *str;
- {
- register struct termentry *t = term;
-
- if (clip_point(x, y)) return;
-
- (*t->put_text)(x,y,str);
- }
-
- /* seems sensible to put the justification in here too..? */
- void clip_put_text_just(x,y,str,just)
- unsigned int x,y;
- char *str;
- int just;
- {
- register struct termentry *t = term;
- if (clip_point(x,y)) return;
- if (!(*t->justify_text)(just)) {
- assert(CENTRE==1 && RIGHT==2);
- x -= (t->h_char * strlen(str) * just)/2;
- }
- (*t->put_text)(x,y,str);
- }
-
- /* single edge intersection algorithm */
- /* Given two points, one inside and one outside the plot, return
- * the point where an edge of the plot intersects the line segment defined
- * by the two points.
- */
- void edge3d_intersect(points, i, ex, ey, ez)
- struct coordinate GPHUGE *points; /* the points array */
- int i; /* line segment from point i-1 to point i */
- double *ex, *ey, *ez; /* the point where it crosses an edge */
- {
- /* global x_min3d, x_max3d, y_min3d, y_max3d, min3d_z, max3d_z */
- int count;
- double ix = points[i-1].x;
- double iy = points[i-1].y;
- double iz = points[i-1].z;
- double ox = points[i].x;
- double oy = points[i].y;
- double oz = points[i].z;
- double x, y, z; /* possible intersection point */
-
- if(points[i].type == INRANGE)
- {
- /* swap points around so that ix/ix/iz are INRANGE and ox/oy/oz are OUTRANGE */
- x = ix;ix = ox;ox = x;
- y = iy;iy = oy;oy = y;
- z = iz;iz = oz;oz = z;
- }
-
- /* nasty degenerate cases, effectively drawing to an infinity point (?)
- cope with them here, so don't process them as a "real" OUTRANGE point
-
- If more than one coord is -VERYLARGE, then can't ratio the "infinities"
- so drop out by returning the INRANGE point.
-
- Obviously, only need to test the OUTRANGE point (coordinates) */
-
- /* nasty degenerate cases, effectively drawing to an infinity point (?)
- cope with them here, so don't process them as a "real" OUTRANGE point
-
- If more than one coord is -VERYLARGE, then can't ratio the "infinities"
- so drop out by returning FALSE */
-
- count = 0;
- if(ox == -VERYLARGE) count++;
- if(oy == -VERYLARGE) count++;
- if(oz == -VERYLARGE) count++;
-
- /* either doesn't pass through 3D volume *or*
- can't ratio infinities to get a direction to draw line, so return the INRANGE point */
- if(count > 1){
- *ex = ix;
- *ey = iy;
- *ez = iz;
-
- return;
- }
-
- if(count == 1)
- {
- *ex = ix;
- *ey = iy;
- *ez = iz;
-
- if(ox == -VERYLARGE)
- {
- *ex = x_min3d;
- return;
- }
-
- if(oy == -VERYLARGE)
- {
- *ey = y_min3d;
- return;
- }
-
- /* obviously oz is -VERYLARGE and (ox != -VERYLARGE && oy != -VERYLARGE) */
- *ez = min3d_z;
- return;
- }
-
- /*
- * Can't have case (ix == ox && iy == oy && iz == oz) as one point
- * is INRANGE and one point is OUTRANGE.
- */
- if(ix == ox) {
- if(iy == oy) {
- /* line parallel to z axis */
-
- /* assume inrange(iy, y_min3d, y_max3d) && inrange(ix, x_min3d, x_max3d) */
- *ex = ix; /* == ox */
- *ey = iy; /* == oy */
-
- if (inrange(max3d_z, iz, oz))
- *ez = max3d_z;
- else if (inrange(min3d_z, iz, oz))
- *ez = min3d_z;
- else {
- graph_error("error in edge3d_intersect");
- }
-
- return;
- }
-
- if(iz == oz) {
- /* line parallel to y axis */
-
- /* assume inrange(iz, min3d_z, max3d_z) && inrange(ix, x_min3d, x_max3d) */
- *ex = ix; /* == ox */
- *ez = iz; /* == oz */
-
- if (inrange(y_max3d, iy, oy))
- *ey = y_max3d;
- else if (inrange(y_min3d, iy, oy))
- *ey = y_min3d;
- else {
- graph_error("error in edge3d_intersect");
- }
-
- return;
- }
-
- /* nasty 2D slanted line in a yz plane */
-
- /* does it intersect y_min3d edge */
- if (inrange(y_min3d, iy, oy) && y_min3d != iy && y_min3d != oy) {
- z = iz + (y_min3d-iy) * ((oz-iz) / (oy-iy));
- if (inrange(z, min3d_z, max3d_z)) {
- *ex = ix;
- *ey = y_min3d;
- *ez = z;
- return;
- }
- }
-
- /* does it intersect y_max3d edge */
- if (inrange(y_max3d, iy, oy) && y_max3d != iy && y_max3d != oy) {
- z = iz + (y_max3d-iy) * ((oz-iz) / (oy-iy));
- if (inrange(z, min3d_z, max3d_z)) {
- *ex = ix;
- *ey = y_max3d;
- *ez = z;
- return;
- }
- }
-
- /* does it intersect min3d_z edge */
- if (inrange(min3d_z, iz, oz) && min3d_z != iz && min3d_z != oz) {
- y = iy + (min3d_z-iz) * ((oy-iy) / (oz-iz));
- if (inrange(y, y_min3d, y_max3d)) {
- *ex = ix;
- *ey = y;
- *ez = min3d_z;
- return;
- }
- }
-
- /* does it intersect max3d_z edge */
- if (inrange(max3d_z, iz, oz) && max3d_z != iz && max3d_z != oz) {
- y = iy + (max3d_z-iz) * ((oy-iy) / (oz-iz));
- if (inrange(y, y_min3d, y_max3d)) {
- *ex = ix;
- *ey = y;
- *ez = max3d_z;
- return;
- }
- }
- }
-
- if(iy == oy) {
- /* already checked case (ix == ox && iy == oy) */
- if(oz == iz) {
- /* line parallel to x axis */
-
- /* assume inrange(iz, min3d_z, max3d_z) && inrange(iy, y_min3d, y_max3d) */
- *ey = iy; /* == oy */
- *ez = iz; /* == oz */
-
- if (inrange(x_max3d, ix, ox))
- *ex = x_max3d;
- else if (inrange(x_min3d, ix, ox))
- *ex = x_min3d;
- else {
- graph_error("error in edge3d_intersect");
- }
-
- return;
- }
-
- /* nasty 2D slanted line in an xz plane */
-
- /* does it intersect x_min3d edge */
- if (inrange(x_min3d, ix, ox) && x_min3d != ix && x_min3d != ox) {
- z = iz + (x_min3d-ix) * ((oz-iz) / (ox-ix));
- if (inrange(z, min3d_z, max3d_z)) {
- *ex = x_min3d;
- *ey = iy;
- *ez = z;
- return;
- }
- }
-
- /* does it intersect x_max3d edge */
- if (inrange(x_max3d, ix, ox) && x_max3d != ix && x_max3d != ox) {
- z = iz + (x_max3d-ix) * ((oz-iz) / (ox-ix));
- if (inrange(z, min3d_z, max3d_z)) {
- *ex = x_max3d;
- *ey = iy;
- *ez = z;
- return;
- }
- }
-
- /* does it intersect min3d_z edge */
- if (inrange(min3d_z, iz, oz) && min3d_z != iz && min3d_z != oz) {
- x = ix + (min3d_z-iz) * ((ox-ix) / (oz-iz));
- if (inrange(x, x_min3d, x_max3d)) {
- *ex = x;
- *ey = iy;
- *ez = min3d_z;
- return;
- }
- }
-
- /* does it intersect max3d_z edge */
- if (inrange(max3d_z, iz, oz) && max3d_z != iz && max3d_z != oz) {
- x = ix + (max3d_z-iz) * ((ox-ix) / (oz-iz));
- if (inrange(x, x_min3d, x_max3d)) {
- *ex = x;
- *ey = iy;
- *ez = max3d_z;
- return;
- }
- }
- }
-
- if(iz == oz) {
- /* already checked cases (ix == ox && iz == oz) and (iy == oy && iz == oz) */
-
- /* nasty 2D slanted line in an xy plane */
-
- /* assume inrange(oz, min3d_z, max3d_z) */
-
- /* does it intersect x_min3d edge */
- if (inrange(x_min3d, ix, ox) && x_min3d != ix && x_min3d != ox) {
- y = iy + (x_min3d-ix) * ((oy-iy) / (ox-ix));
- if (inrange(y, y_min3d, y_max3d)) {
- *ex = x_min3d;
- *ey = y;
- *ez = iz;
- return;
- }
- }
-
- /* does it intersect x_max3d edge */
- if (inrange(x_max3d, ix, ox) && x_max3d != ix && x_max3d != ox) {
- y = iy + (x_max3d-ix) * ((oy-iy) / (ox-ix));
- if (inrange(y, y_min3d, y_max3d)) {
- *ex = x_max3d;
- *ey = y;
- *ez = iz;
- return;
- }
- }
-
- /* does it intersect y_min3d edge */
- if (inrange(y_min3d, iy, oy) && y_min3d != iy && y_min3d != oy) {
- x = ix + (y_min3d-iy) * ((ox-ix) / (oy-iy));
- if (inrange(x, x_min3d, x_max3d)) {
- *ex = x;
- *ey = y_min3d;
- *ez = iz;
- return;
- }
- }
-
- /* does it intersect y_max3d edge */
- if (inrange(y_max3d, iy, oy) && y_max3d != iy && y_max3d != oy) {
- x = ix + (y_max3d-iy) * ((ox-ix) / (oy-iy));
- if (inrange(x, x_min3d, x_max3d)) {
- *ex = x;
- *ey = y_max3d;
- *ez = iz;
- return;
- }
- }
- }
-
- /* really nasty general slanted 3D case */
-
- /* does it intersect x_min3d edge */
- if (inrange(x_min3d, ix, ox) && x_min3d != ix && x_min3d != ox) {
- y = iy + (x_min3d-ix) * ((oy-iy) / (ox-ix));
- z = iz + (x_min3d-ix) * ((oz-iz) / (ox-ix));
- if (inrange(y, y_min3d, y_max3d) && inrange(z, min3d_z, max3d_z)) {
- *ex = x_min3d;
- *ey = y;
- *ez = z;
- return;
- }
- }
-
- /* does it intersect x_max3d edge */
- if (inrange(x_max3d, ix, ox) && x_max3d != ix && x_max3d != ox) {
- y = iy + (x_max3d-ix) * ((oy-iy) / (ox-ix));
- z = iz + (x_max3d-ix) * ((oz-iz) / (ox-ix));
- if (inrange(y, y_min3d, y_max3d) && inrange(z, min3d_z, max3d_z)) {
- *ex = x_max3d;
- *ey = y;
- *ez = z;
- return;
- }
- }
-
- /* does it intersect y_min3d edge */
- if (inrange(y_min3d, iy, oy) && y_min3d != iy && y_min3d != oy) {
- x = ix + (y_min3d-iy) * ((ox-ix) / (oy-iy));
- z = iz + (y_min3d-iy) * ((oz-iz) / (oy-iy));
- if (inrange(x, x_min3d, x_max3d) && inrange(z, min3d_z, max3d_z)) {
- *ex = x;
- *ey = y_min3d;
- *ez = z;
- return;
- }
- }
-
- /* does it intersect y_max3d edge */
- if (inrange(y_max3d, iy, oy) && y_max3d != iy && y_max3d != oy) {
- x = ix + (y_max3d-iy) * ((ox-ix) / (oy-iy));
- z = iz + (y_max3d-iy) * ((oz-iz) / (oy-iy));
- if (inrange(x, x_min3d, x_max3d) && inrange(z, min3d_z, max3d_z)) {
- *ex = x;
- *ey = y_max3d;
- *ez = z;
- return;
- }
- }
-
- /* does it intersect min3d_z edge */
- if (inrange(min3d_z, iz, oz) && min3d_z != iz && min3d_z != oz) {
- x = ix + (min3d_z-iz) * ((ox-ix) / (oz-iz));
- y = iy + (min3d_z-iz) * ((oy-iy) / (oz-iz));
- if (inrange(x, x_min3d, x_max3d) && inrange(y, y_min3d, y_max3d)) {
- *ex = x;
- *ey = y;
- *ez = min3d_z;
- return;
- }
- }
-
- /* does it intersect max3d_z edge */
- if (inrange(max3d_z, iz, oz) && max3d_z != iz && max3d_z != oz) {
- x = ix + (max3d_z-iz) * ((ox-ix) / (oz-iz));
- y = iy + (max3d_z-iz) * ((oy-iy) / (oz-iz));
- if (inrange(x, x_min3d, x_max3d) && inrange(y, y_min3d, y_max3d)) {
- *ex = x;
- *ey = y;
- *ez = max3d_z;
- return;
- }
- }
-
- /* If we reach here, the inrange point is on the edge, and
- * the line segment from the outrange point does not cross any
- * other edges to get there. In this case, we return the inrange
- * point as the 'edge' intersection point. This will basically draw
- * line.
- */
- *ex = ix;
- *ey = iy;
- *ez = iz;
- return;
- }
-
- /* double edge intersection algorithm */
- /* Given two points, both outside the plot, return
- * the points where an edge of the plot intersects the line segment defined
- * by the two points. There may be zero, one, two, or an infinite number
- * of intersection points. (One means an intersection at a corner, infinite
- * means overlaying the edge itself). We return FALSE when there is nothing
- * to draw (zero intersections), and TRUE when there is something to
- * draw (the one-point case is a degenerate of the two-point case and we do
- * not distinguish it - we draw it anyway).
- */
- TBOOLEAN /* any intersection? */
- two_edge3d_intersect(points, i, lx, ly, lz)
- struct coordinate GPHUGE *points; /* the points array */
- int i; /* line segment from point i-1 to point i */
- double *lx, *ly, *lz; /* lx[2], ly[2], lz[2]: points where it crosses edges */
- {
- int count;
- /* global x_min3d, x_max3d, y_min3d, y_max3d, min3d_z, max3d_z */
- double ix = points[i-1].x;
- double iy = points[i-1].y;
- double iz = points[i-1].z;
- double ox = points[i].x;
- double oy = points[i].y;
- double oz = points[i].z;
- double t[6];
- double swap;
- double x, y, z; /* possible intersection point */
- double t_min, t_max;
-
- /* nasty degenerate cases, effectively drawing to an infinity point (?)
- cope with them here, so don't process them as a "real" OUTRANGE point
-
- If more than one coord is -VERYLARGE, then can't ratio the "infinities"
- so drop out by returning FALSE */
-
- count = 0;
- if(ix == -VERYLARGE) count++;
- if(ox == -VERYLARGE) count++;
- if(iy == -VERYLARGE) count++;
- if(oy == -VERYLARGE) count++;
- if(iz == -VERYLARGE) count++;
- if(oz == -VERYLARGE) count++;
-
- /* either doesn't pass through 3D volume *or*
- can't ratio infinities to get a direction to draw line, so simply return(FALSE) */
- if(count > 1){
- return(FALSE);
- }
-
- if(ox == -VERYLARGE || ix == -VERYLARGE)
- {
- if(ix == -VERYLARGE)
- {
- /* swap points so ix/iy/iz don't have a -VERYLARGE component */
- x = ix;ix = ox;ox = x;
- y = iy;iy = oy;oy = y;
- z = iz;iz = oz;oz = z;
- }
-
- /* check actually passes through the 3D graph volume */
- if(ix > x_max3d && inrange(iy, y_min3d, y_max3d) && inrange(iz, min3d_z, max3d_z))
- {
- lx[0] = x_min3d;
- ly[0] = iy;
- lz[0] = iz;
-
- lx[1] = x_max3d;
- ly[1] = iy;
- lz[1] = iz;
-
- return(TRUE);
- }
- else {
- return(FALSE);
- }
- }
-
- if(oy == -VERYLARGE || iy == -VERYLARGE)
- {
- if(iy == -VERYLARGE)
- {
- /* swap points so ix/iy/iz don't have a -VERYLARGE component */
- x = ix;ix = ox;ox = x;
- y = iy;iy = oy;oy = y;
- z = iz;iz = oz;oz = z;
- }
-
- /* check actually passes through the 3D graph volume */
- if(iy > y_max3d && inrange(ix, x_min3d, x_max3d) && inrange(iz, min3d_z, max3d_z))
- {
- lx[0] = ix;
- ly[0] = y_min3d;
- lz[0] = iz;
-
- lx[1] = ix;
- ly[1] = y_max3d;
- lz[1] = iz;
-
- return(TRUE);
- }
- else {
- return(FALSE);
- }
- }
-
- if(oz == -VERYLARGE || iz == -VERYLARGE)
- {
- if(iz == -VERYLARGE)
- {
- /* swap points so ix/iy/iz don't have a -VERYLARGE component */
- x = ix;ix = ox;ox = x;
- y = iy;iy = oy;oy = y;
- z = iz;iz = oz;oz = z;
- }
-
- /* check actually passes through the 3D graph volume */
- if(iz > max3d_z && inrange(ix, x_min3d, x_max3d) && inrange(iy, y_min3d, y_max3d))
- {
- lx[0] = ix;
- ly[0] = iy;
- lz[0] = min3d_z;
-
- lx[1] = ix;
- ly[1] = iy;
- lz[1] = max3d_z;
-
- return(TRUE);
- }
- else {
- return(FALSE);
- }
- }
-
- /*
- * Quick outcode tests on the 3d graph volume
- */
-
- /*
- * test z coord first --- most surface OUTRANGE points generated between
- * min3d_z and z_min3d (i.e. when ticslevel is non-zero)
- */
- if(GPMAX(iz,oz) < min3d_z || GPMIN(iz,oz) > max3d_z)
- return(FALSE);
-
- if(GPMAX(ix,ox) < x_min3d || GPMIN(ix,ox) > x_max3d)
- return(FALSE);
-
- if(GPMAX(iy,oy) < y_min3d || GPMIN(iy,oy) > y_max3d)
- return(FALSE);
-
- /*
- * Special horizontal/vertical, etc. cases are checked and remaining
- * slant lines are checked separately.
- *
- * The slant line intersections are solved using the parametric form
- * of the equation for a line, since if we test x/y/z min/max planes explicitly
- * then e.g. a line passing through a corner point (x_min,y_min,z_min)
- * actually intersects all 3 planes and hence further tests would be required
- * to anticipate this and similar situations.
- */
-
- /*
- * Can have case (ix == ox && iy == oy && iz == oz) as both points OUTRANGE
- */
- if(ix == ox && iy == oy && iz == oz)
- {
- /* but as only define single outrange point, can't intersect 3D graph volume */
- return(FALSE);
- }
-
- if(ix == ox) {
- if(iy == oy) {
- /* line parallel to z axis */
-
- /* x and y coords must be in range, and line must span both min3d_z and max3d_z */
- /* note that spanning min3d_z implies spanning max3d_z as both points OUTRANGE */
- if(!inrange(ix, x_min3d, x_max3d) || !inrange(iy, y_min3d, y_max3d))
- {
- return(FALSE);
- }
-
- if (inrange(min3d_z, iz, oz)) {
- lx[0] = ix;
- ly[0] = iy;
- lz[0] = min3d_z;
-
- lx[1] = ix;
- ly[1] = iy;
- lz[1] = max3d_z;
-
- return(TRUE);
- } else
- return(FALSE);
- }
-
- if(iz == oz) {
- /* line parallel to y axis */
-
- /* x and z coords must be in range, and line must span both y_min3d and y_max3d */
- /* note that spanning y_min3d implies spanning y_max3d, as both points OUTRANGE */
- if(!inrange(ix, x_min3d, x_max3d) || !inrange(iz, min3d_z, max3d_z))
- {
- return(FALSE);
- }
-
- if (inrange(y_min3d, iy, oy)) {
- lx[0] = ix;
- ly[0] = y_min3d;
- lz[0] = iz;
-
- lx[1] = ix;
- ly[1] = y_max3d;
- lz[1] = iz;
-
- return(TRUE);
- } else
- return(FALSE);
- }
-
- /* nasty 2D slanted line in a yz plane */
-
- if(!inrange(ox, x_min3d, x_max3d))
- return(FALSE);
-
- t[0] = (y_min3d - iy)/(oy - iy);
- t[1] = (y_max3d - iy)/(oy - iy);
-
- if(t[0] > t[1]) {
- swap = t[0];t[0] = t[1];t[1] = swap;
- }
-
- t[2] = (min3d_z - iz)/(oz - iz);
- t[3] = (max3d_z - iz)/(oz - iz);
-
- if(t[2] > t[3]) {
- swap = t[2];t[2] = t[3];t[3] = swap;
- }
-
- t_min = GPMAX(GPMAX(t[0],t[2]),0.0);
- t_max = GPMIN(GPMIN(t[1],t[3]),1.0);
-
- if(t_min > t_max)
- return(FALSE);
-
- lx[0] = ix;
- ly[0] = iy + t_min * (oy - iy);
- lz[0] = iz + t_min * (oz - iz);
-
- lx[1] = ix;
- ly[1] = iy + t_max * (oy - iy);
- lz[1] = iz + t_max * (oz - iz);
-
- /*
- * Can only have 0 or 2 intersection points -- only need test one coord
- */
- if(inrange(ly[0], y_min3d, y_max3d) &&
- inrange(lz[0], min3d_z, max3d_z))
- {
- return(TRUE);
- }
-
- return(FALSE);
- }
-
- if(iy == oy) {
- /* already checked case (ix == ox && iy == oy) */
- if(oz == iz) {
- /* line parallel to x axis */
-
- /* y and z coords must be in range, and line must span both x_min3d and x_max3d */
- /* note that spanning x_min3d implies spanning x_max3d, as both points OUTRANGE */
- if(!inrange(iy, y_min3d, y_max3d) || !inrange(iz, min3d_z, max3d_z))
- {
- return(FALSE);
- }
-
- if (inrange(x_min3d, ix, ox)) {
- lx[0] = x_min3d;
- ly[0] = iy;
- lz[0] = iz;
-
- lx[1] = x_max3d;
- ly[1] = iy;
- lz[1] = iz;
-
- return(TRUE);
- } else
- return(FALSE);
- }
-
- /* nasty 2D slanted line in an xz plane */
-
- if(!inrange(oy, y_min3d, y_max3d))
- return(FALSE);
-
- t[0] = (x_min3d - ix)/(ox - ix);
- t[1] = (x_max3d - ix)/(ox - ix);
-
- if(t[0] > t[1]) {
- swap = t[0];t[0] = t[1];t[1] = swap;
- }
-
- t[2] = (min3d_z - iz)/(oz - iz);
- t[3] = (max3d_z - iz)/(oz - iz);
-
- if(t[2] > t[3]) {
- swap = t[2];t[2] = t[3];t[3] = swap;
- }
-
- t_min = GPMAX(GPMAX(t[0],t[2]),0.0);
- t_max = GPMIN(GPMIN(t[1],t[3]),1.0);
-
- if(t_min > t_max)
- return(FALSE);
-
- lx[0] = ix + t_min * (ox - ix);
- ly[0] = iy;
- lz[0] = iz + t_min * (oz - iz);
-
- lx[1] = ix + t_max * (ox - ix);
- ly[1] = iy;
- lz[1] = iz + t_max * (oz - iz);
-
- /*
- * Can only have 0 or 2 intersection points -- only need test one coord
- */
- if(inrange(lx[0], x_min3d, x_max3d) &&
- inrange(lz[0], min3d_z, max3d_z))
- {
- return(TRUE);
- }
-
- return(FALSE);
- }
-
- if(iz == oz) {
- /* already checked cases (ix == ox && iz == oz) and (iy == oy && iz == oz) */
-
- /* nasty 2D slanted line in an xy plane */
-
- if(!inrange(oz, min3d_z, max3d_z))
- return(FALSE);
-
- t[0] = (x_min3d - ix)/(ox - ix);
- t[1] = (x_max3d - ix)/(ox - ix);
-
- if(t[0] > t[1]) {
- swap = t[0];t[0] = t[1];t[1] = swap;
- }
-
- t[2] = (y_min3d - iy)/(oy - iy);
- t[3] = (y_max3d - iy)/(oy - iy);
-
- if(t[2] > t[3]) {
- swap = t[2];t[2] = t[3];t[3] = swap;
- }
-
- t_min = GPMAX(GPMAX(t[0],t[2]),0.0);
- t_max = GPMIN(GPMIN(t[1],t[3]),1.0);
-
- if(t_min > t_max)
- return(FALSE);
-
- lx[0] = ix + t_min * (ox - ix);
- ly[0] = iy + t_min * (oy - iy);
- lz[0] = iz;
-
- lx[1] = ix + t_max * (ox - ix);
- ly[1] = iy + t_max * (oy - iy);
- lz[1] = iz;
-
- /*
- * Can only have 0 or 2 intersection points -- only need test one coord
- */
- if(inrange(lx[0], x_min3d, x_max3d) &&
- inrange(ly[0], y_min3d, y_max3d))
- {
- return(TRUE);
- }
-
- return(FALSE);
- }
-
- /* really nasty general slanted 3D case */
-
- /*
- Solve parametric equation
-
- (ix, iy, iz) + t (diff_x, diff_y, diff_z)
-
- where 0.0 <= t <= 1.0 and
-
- diff_x = (ox - ix);
- diff_y = (oy - iy);
- diff_z = (oz - iz);
- */
-
- t[0] = (x_min3d - ix)/(ox - ix);
- t[1] = (x_max3d - ix)/(ox - ix);
-
- if(t[0] > t[1]) {
- swap = t[0];t[0] = t[1];t[1] = swap;
- }
-
- t[2] = (y_min3d - iy)/(oy - iy);
- t[3] = (y_max3d - iy)/(oy - iy);
-
- if(t[2] > t[3]) {
- swap = t[2];t[2] = t[3];t[3] = swap;
- }
-
- t[4] = (iz == oz) ? 0.0 : (min3d_z - iz)/(oz - iz);
- t[5] = (iz == oz) ? 1.0 : (max3d_z - iz)/(oz - iz);
-
- if(t[4] > t[5]) {
- swap = t[4];t[4] = t[5];t[5] = swap;
- }
-
- t_min = GPMAX(GPMAX(t[0],t[2]),GPMAX(t[4],0.0));
- t_max = GPMIN(GPMIN(t[1],t[3]),GPMIN(t[5],1.0));
-
- if(t_min > t_max)
- return(FALSE);
-
- lx[0] = ix + t_min * (ox - ix);
- ly[0] = iy + t_min * (oy - iy);
- lz[0] = iz + t_min * (oz - iz);
-
- lx[1] = ix + t_max * (ox - ix);
- ly[1] = iy + t_max * (oy - iy);
- lz[1] = iz + t_max * (oz - iz);
-
- /*
- * Can only have 0 or 2 intersection points -- only need test one coord
- */
- if(inrange(lx[0], x_min3d, x_max3d) &&
- inrange(ly[0], y_min3d, y_max3d) &&
- inrange(lz[0], min3d_z, max3d_z))
- {
- return(TRUE);
- }
-
- return(FALSE);
- }
-
-
-
-