git.maemo.org Git - gnuplot/blob - src/util3d.c

   1 #ifndef lint
   2 static char *RCSid() { return RCSid("$Id: util3d.c,v 1.27.2.4 2009/01/07 22:58:27 sfeam Exp $"); }
   3 #endif
   4
   5 /* GNUPLOT - util3d.c */
   6
   7 /*[
   8  * Copyright 1986 - 1993, 1998, 2004   Thomas Williams, Colin Kelley
   9  *
  10  * Permission to use, copy, and distribute this software and its
  11  * documentation for any purpose with or without fee is hereby granted,
  12  * provided that the above copyright notice appear in all copies and
  13  * that both that copyright notice and this permission notice appear
  14  * in supporting documentation.
  15  *
  16  * Permission to modify the software is granted, but not the right to
  17  * distribute the complete modified source code.  Modifications are to
  18  * be distributed as patches to the released version.  Permission to
  19  * distribute binaries produced by compiling modified sources is granted,
  20  * provided you
  21  *   1. distribute the corresponding source modifications from the
  22  *    released version in the form of a patch file along with the binaries,
  23  *   2. add special version identification to distinguish your version
  24  *    in addition to the base release version number,
  25  *   3. provide your name and address as the primary contact for the
  26  *    support of your modified version, and
  27  *   4. retain our contact information in regard to use of the base
  28  *    software.
  29  * Permission to distribute the released version of the source code along
  30  * with corresponding source modifications in the form of a patch file is
  31  * granted with same provisions 2 through 4 for binary distributions.
  32  *
  33  * This software is provided "as is" without express or implied warranty
  34  * to the extent permitted by applicable law.
  35 ]*/
  36
  37
  38 /*
  39  * 19 September 1992  Lawrence Crowl  (crowl@cs.orst.edu)
  40  * Added user-specified bases for log scaling.
  41  *
  42  * 3.6 - split graph3d.c into graph3d.c (graph),
  43  *                            util3d.c (intersections, etc)
  44  *                            hidden3d.c (hidden-line removal code)
  45  *
  46  */
  47
  48 #include "util3d.h"
  49
  50 #include "axis.h"
  51 #include "hidden3d.h"
  52 #include "pm3d.h"
  53 #include "term_api.h"
  54
  55 /* HBB 990826: all that stuff referenced from other modules is now
  56  * exported in graph3d.h, instead of being listed here */
  57
  58 /* Prototypes for local functions */
  59 static void mat_unit __PROTO((transform_matrix mat));
  60 static GP_INLINE void draw3d_point_unconditional __PROTO((p_vertex, struct lp_style_type *));
  61
  62 static void
  63 mat_unit(transform_matrix mat)
  64 {
  65     int i, j;
  66
  67     for (i = 0; i < 4; i++)
  68         for (j = 0; j < 4; j++)
  69             if (i == j)
  70                 mat[i][j] = 1.0;
  71             else
  72                 mat[i][j] = 0.0;
  73 }
  74
  75 #if 0 /* HBB 990829: unused --> commented out */
  76 void
  77 mat_trans(double tx, double ty, double tz, transform_matrix mat)
  78 {
  79     mat_unit(mat);              /* Make it unit matrix. */
  80     mat[3][0] = tx;
  81     mat[3][1] = ty;
  82     mat[3][2] = tz;
  83 }
  84 #endif /* commented out */
  85
  86 void
  87 mat_scale(double sx, double sy, double sz, transform_matrix mat)
  88 {
  89     mat_unit(mat);              /* Make it unit matrix. */
  90     mat[0][0] = sx;
  91     mat[1][1] = sy;
  92     mat[2][2] = sz;
  93 }
  94
  95 void
  96 mat_rot_x(double teta, transform_matrix mat)
  97 {
  98     double cos_teta, sin_teta;
  99
 100     teta *= DEG2RAD;
 101     cos_teta = cos(teta);
 102     sin_teta = sin(teta);
 103
 104     mat_unit(mat);              /* Make it unit matrix. */
 105     mat[1][1] = cos_teta;
 106     mat[1][2] = -sin_teta;
 107     mat[2][1] = sin_teta;
 108     mat[2][2] = cos_teta;
 109 }
 110
 111 #if 0 /* HBB 990829: unused --> commented out */
 112 void
 113 mat_rot_y(double teta, transform_matrix mat)
 114 {
 115     double cos_teta, sin_teta;
 116
 117     teta *= DEG2RAD;
 118     cos_teta = cos(teta);
 119     sin_teta = sin(teta);
 120
 121     mat_unit(mat);              /* Make it unit matrix. */
 122     mat[0][0] = cos_teta;
 123     mat[0][2] = -sin_teta;
 124     mat[2][0] = sin_teta;
 125     mat[2][2] = cos_teta;
 126 }
 127 #endif /* commented out */
 128
 129 void
 130 mat_rot_z(double teta, transform_matrix mat)
 131 {
 132     double cos_teta, sin_teta;
 133
 134     teta *= DEG2RAD;
 135     cos_teta = cos(teta);
 136     sin_teta = sin(teta);
 137
 138     mat_unit(mat);              /* Make it unit matrix. */
 139     mat[0][0] = cos_teta;
 140     mat[0][1] = -sin_teta;
 141     mat[1][0] = sin_teta;
 142     mat[1][1] = cos_teta;
 143 }
 144
 145 /* Multiply two transform_matrix. Result can be one of two operands. */
 146 void
 147 mat_mult(
 148     transform_matrix mat_res,
 149     transform_matrix mat1, transform_matrix mat2)
 150 {
 151     int i, j, k;
 152     transform_matrix mat_res_temp;
 153
 154     for (i = 0; i < 4; i++)
 155         for (j = 0; j < 4; j++) {
 156             mat_res_temp[i][j] = 0;
 157             for (k = 0; k < 4; k++)
 158                 mat_res_temp[i][j] += mat1[i][k] * mat2[k][j];
 159         }
 160     for (i = 0; i < 4; i++)
 161         for (j = 0; j < 4; j++)
 162             mat_res[i][j] = mat_res_temp[i][j];
 163 }
 164
 165 #define IN_AXIS_RANGE(val, axis)                                        \
 166         inrange((val), axis_array[axis].min, axis_array[axis].max)
 167
 168
 169 /* single edge intersection algorithm */
 170 /* Given two points, one inside and one outside the plot, return
 171  * the point where an edge of the plot intersects the line segment defined
 172  * by the two points.
 173  */
 174 void
 175 edge3d_intersect(
 176     struct coordinate GPHUGE *points,   /* the points array */
 177     int i,                              /* line segment from point i-1 to point i */
 178     double *ex, double *ey, double *ez) /* the point where it crosses an edge */
 179 {
 180     int count;
 181     double ix = points[i - 1].x;
 182     double iy = points[i - 1].y;
 183     double iz = points[i - 1].z;
 184     double ox = points[i].x;
 185     double oy = points[i].y;
 186     double oz = points[i].z;
 187     double x, y, z;             /* possible intersection point */
 188
 189     if (points[i].type == INRANGE) {
 190         /* swap points around so that ix/ix/iz are INRANGE and ox/oy/oz are OUTRANGE */
 191         x = ix;
 192         ix = ox;
 193         ox = x;
 194         y = iy;
 195         iy = oy;
 196         oy = y;
 197         z = iz;
 198         iz = oz;
 199         oz = z;
 200     }
 201
 202     /* nasty degenerate cases, effectively drawing to an infinity point (?)
 203        cope with them here, so don't process them as a "real" OUTRANGE point
 204
 205        If more than one coord is -VERYLARGE, then can't ratio the "infinities"
 206        so drop out by returning FALSE */
 207
 208     count = 0;
 209     if (ox == -VERYLARGE)
 210         count++;
 211     if (oy == -VERYLARGE)
 212         count++;
 213     if (oz == -VERYLARGE)
 214         count++;
 215
 216     /* either doesn't pass through 3D volume *or*
 217        can't ratio infinities to get a direction to draw line, so return the INRANGE point */
 218     if (count > 1) {
 219         *ex = ix;
 220         *ey = iy;
 221         *ez = iz;
 222
 223         return;
 224     }
 225     if (count == 1) {
 226         *ex = ix;
 227         *ey = iy;
 228         *ez = iz;
 229
 230         if (ox == -VERYLARGE) {
 231             *ex = AXIS_ACTUAL_MIN(FIRST_X_AXIS);
 232             return;
 233         }
 234         if (oy == -VERYLARGE) {
 235             *ey = AXIS_ACTUAL_MIN(FIRST_Y_AXIS);
 236             return;
 237         }
 238         /* obviously oz is -VERYLARGE and (ox != -VERYLARGE && oy != -VERYLARGE) */
 239         *ez = AXIS_ACTUAL_MIN(FIRST_Z_AXIS);
 240         return;
 241     }
 242     /*
 243      * Can't have case (ix == ox && iy == oy && iz == oz) as one point
 244      * is INRANGE and one point is OUTRANGE.
 245      */
 246     if (ix == ox) {
 247         if (iy == oy) {
 248             /* line parallel to z axis */
 249
 250             /* assume iy in yrange, && ix in xrange */
 251             *ex = ix;           /* == ox */
 252             *ey = iy;           /* == oy */
 253
 254             if (inrange(AXIS_ACTUAL_MAX(FIRST_Z_AXIS), iz, oz))
 255                 *ez = AXIS_ACTUAL_MAX(FIRST_Z_AXIS);
 256             else if (inrange(AXIS_ACTUAL_MIN(FIRST_Z_AXIS), iz, oz))
 257                 *ez = AXIS_ACTUAL_MIN(FIRST_Z_AXIS);
 258             else {
 259                 graph_error("error in edge3d_intersect");
 260             }
 261
 262             return;
 263         }
 264         if (iz == oz) {
 265             /* line parallel to y axis */
 266
 267             /* assume iz in zrange && ix in xrange */
 268             *ex = ix;           /* == ox */
 269             *ez = iz;           /* == oz */
 270
 271             if (inrange(AXIS_ACTUAL_MAX(FIRST_Y_AXIS), iy, oy))
 272                 *ey = AXIS_ACTUAL_MAX(FIRST_Y_AXIS);
 273             else if (inrange(AXIS_ACTUAL_MIN(FIRST_Y_AXIS), iy, oy))
 274                 *ey = AXIS_ACTUAL_MIN(FIRST_Y_AXIS);
 275             else {
 276                 graph_error("error in edge3d_intersect");
 277             }
 278
 279             return;
 280         }
 281
 282         /* nasty 2D slanted line in a yz plane */
 283
 284
 285 #define INTERSECT_PLANE(cut, axis, eff, eff_axis, res_x, res_y, res_z)  \
 286         do {                                                            \
 287             if (inrange(cut, i##axis, o##axis)                          \
 288                 && cut != i##axis                                       \
 289                 && cut != o##axis) {                                    \
 290                 eff = (cut - i##axis)                                   \
 291                     * ((o##eff - i##eff) / (o##axis - i##axis))         \
 292                     + i##eff;                                           \
 293                 if (IN_AXIS_RANGE(eff, eff_axis)) {                     \
 294                     *ex = res_x;                                        \
 295                     *ey = res_y;                                        \
 296                     *ez = res_z;                                        \
 297                     return;                                             \
 298                 }                                                       \
 299             }                                                           \
 300         } while (0)
 301
 302         INTERSECT_PLANE(AXIS_ACTUAL_MIN(FIRST_Y_AXIS), y, z, FIRST_Z_AXIS,
 303                         ix, AXIS_ACTUAL_MIN(FIRST_Y_AXIS), z);
 304         INTERSECT_PLANE(AXIS_ACTUAL_MAX(FIRST_Y_AXIS), y, z, FIRST_Z_AXIS,
 305                         ix, AXIS_ACTUAL_MAX(FIRST_Y_AXIS), z);
 306         INTERSECT_PLANE(AXIS_ACTUAL_MIN(FIRST_Z_AXIS), z, y, FIRST_Y_AXIS,
 307                         ix, y, AXIS_ACTUAL_MIN(FIRST_Z_AXIS));
 308         INTERSECT_PLANE(AXIS_ACTUAL_MAX(FIRST_Z_AXIS), z, y, FIRST_Y_AXIS,
 309                         ix, y, AXIS_ACTUAL_MAX(FIRST_Z_AXIS));
 310     } /* if (ix == ox) */
 311
 312     if (iy == oy) {
 313         /* already checked case (ix == ox && iy == oy) */
 314         if (oz == iz) {
 315             /* line parallel to x axis */
 316
 317             /* assume inrange(iz) && inrange(iy) */
 318             *ey = iy;           /* == oy */
 319             *ez = iz;           /* == oz */
 320
 321             if (inrange(AXIS_ACTUAL_MAX(FIRST_X_AXIS), ix, ox))
 322                 *ex = AXIS_ACTUAL_MAX(FIRST_X_AXIS);
 323             else if (inrange(AXIS_ACTUAL_MIN(FIRST_X_AXIS), ix, ox))
 324                 *ex = AXIS_ACTUAL_MIN(FIRST_X_AXIS);
 325             else {
 326                 graph_error("error in edge3d_intersect");
 327             }
 328
 329             return;
 330         }
 331         /* nasty 2D slanted line in an xz plane */
 332
 333         INTERSECT_PLANE(AXIS_ACTUAL_MIN(FIRST_X_AXIS), x, z, FIRST_Z_AXIS,
 334                         AXIS_ACTUAL_MIN(FIRST_X_AXIS), iy, z);
 335         INTERSECT_PLANE(AXIS_ACTUAL_MAX(FIRST_X_AXIS), x, z, FIRST_Z_AXIS,
 336                         AXIS_ACTUAL_MAX(FIRST_X_AXIS), iy, z);
 337         INTERSECT_PLANE(AXIS_ACTUAL_MIN(FIRST_Z_AXIS), z, x, FIRST_X_AXIS,
 338                         x, iy, AXIS_ACTUAL_MIN(FIRST_Z_AXIS));
 339         INTERSECT_PLANE(AXIS_ACTUAL_MAX(FIRST_Z_AXIS), z, x, FIRST_X_AXIS,
 340                         x, iy, AXIS_ACTUAL_MAX(FIRST_Z_AXIS));
 341     } /* if(iy==oy) */
 342
 343     if (iz == oz) {
 344         /* already checked cases (ix == ox && iz == oz) and (iy == oy
 345            && iz == oz) */
 346
 347         /* 2D slanted line in an xy plane */
 348
 349         /* assume inrange(oz) */
 350
 351         INTERSECT_PLANE(AXIS_ACTUAL_MIN(FIRST_X_AXIS), x, y, FIRST_Y_AXIS,
 352                         AXIS_ACTUAL_MIN(FIRST_X_AXIS), y, iz);
 353         INTERSECT_PLANE(AXIS_ACTUAL_MAX(FIRST_X_AXIS), x, y, FIRST_Y_AXIS,
 354                         AXIS_ACTUAL_MAX(FIRST_X_AXIS), y, iz);
 355         INTERSECT_PLANE(AXIS_ACTUAL_MIN(FIRST_Y_AXIS), y, x, FIRST_X_AXIS,
 356                         x, AXIS_ACTUAL_MIN(FIRST_Y_AXIS), iz);
 357         INTERSECT_PLANE(AXIS_ACTUAL_MAX(FIRST_Y_AXIS), y, x, FIRST_X_AXIS,
 358                         x, AXIS_ACTUAL_MAX(FIRST_Y_AXIS), iz);
 359     } /* if(iz==oz) */
 360 #undef INTERSECT_PLANE
 361
 362     /* really nasty general slanted 3D case */
 363
 364 #define INTERSECT_DIAG(cut, axis, eff, eff_axis, eff2, eff2_axis,       \
 365                        res_x, res_y, res_z)                             \
 366         do {                                                            \
 367             if (inrange(cut, i##axis, o##axis)                          \
 368                 && cut != i##axis                                       \
 369                 && cut != o##axis) {                                    \
 370                 eff = (cut - i##axis)                                   \
 371                     * ((o##eff - i##eff) / (o##axis - i##axis))         \
 372                     + i##eff;                                           \
 373                 eff2 = (cut - i##axis)                                  \
 374                     * ((o##eff2 - i##eff2) / (o##axis - i##axis))       \
 375                     + i##eff2;                                          \
 376                 if (IN_AXIS_RANGE(eff, eff_axis)                        \
 377                     && IN_AXIS_RANGE(eff2, eff2_axis)) {                \
 378                     *ex = res_x;                                        \
 379                     *ey = res_y;                                        \
 380                     *ez = res_z;                                        \
 381                     return;                                             \
 382                 }                                                       \
 383             }                                                           \
 384         } while (0)
 385
 386     INTERSECT_DIAG(AXIS_ACTUAL_MIN(FIRST_X_AXIS), x,
 387                    y, FIRST_Y_AXIS, z, FIRST_Z_AXIS,
 388                    AXIS_ACTUAL_MIN(FIRST_X_AXIS), y, z);
 389     INTERSECT_DIAG(AXIS_ACTUAL_MAX(FIRST_X_AXIS), x,
 390                    y, FIRST_Y_AXIS, z, FIRST_Z_AXIS,
 391                    AXIS_ACTUAL_MAX(FIRST_X_AXIS), y, z);
 392
 393     INTERSECT_DIAG(AXIS_ACTUAL_MIN(FIRST_Y_AXIS), y,
 394                    x, FIRST_X_AXIS, z, FIRST_Z_AXIS,
 395                    x, AXIS_ACTUAL_MIN(FIRST_Y_AXIS), z);
 396     INTERSECT_DIAG(AXIS_ACTUAL_MAX(FIRST_Y_AXIS), y,
 397                    x, FIRST_X_AXIS, z, FIRST_Z_AXIS,
 398                    x, AXIS_ACTUAL_MAX(FIRST_Y_AXIS), z);
 399
 400     INTERSECT_DIAG(AXIS_ACTUAL_MIN(FIRST_Z_AXIS), z,
 401                    x, FIRST_X_AXIS, y, FIRST_Y_AXIS,
 402                    x, y, AXIS_ACTUAL_MIN(FIRST_Z_AXIS));
 403     INTERSECT_DIAG(AXIS_ACTUAL_MAX(FIRST_Z_AXIS), z,
 404                    x, FIRST_X_AXIS, y, FIRST_Y_AXIS,
 405                    x, y, AXIS_ACTUAL_MAX(FIRST_Z_AXIS));
 406
 407 #undef INTERSECT_DIAG
 408
 409     /* If we reach here, the inrange point is on the edge, and
 410      * the line segment from the outrange point does not cross any
 411      * other edges to get there. In this case, we return the inrange
 412      * point as the 'edge' intersection point. This will basically draw
 413      * line.
 414      */
 415     *ex = ix;
 416     *ey = iy;
 417     *ez = iz;
 418     return;
 419 }
 420
 421 /* double edge intersection algorithm */
 422 /* Given two points, both outside the plot, return
 423  * the points where an edge of the plot intersects the line segment defined
 424  * by the two points. There may be zero, one, two, or an infinite number
 425  * of intersection points. (One means an intersection at a corner, infinite
 426  * means overlaying the edge itself). We return FALSE when there is nothing
 427  * to draw (zero intersections), and TRUE when there is something to
 428  * draw (the one-point case is a degenerate of the two-point case and we do
 429  * not distinguish it - we draw it anyway).
 430  */
 431 TBOOLEAN                        /* any intersection? */
 432 two_edge3d_intersect(
 433     struct coordinate GPHUGE *points,   /* the points array */
 434     int i,                              /* line segment from point i-1 to point i */
 435     double *lx, double *ly, double *lz) /* lx[2], ly[2], lz[2]: points where it crosses edges */
 436 {
 437     int count;
 438     /* global axis_array[FIRST_{X,Y,Z}_AXIS].{min,max} */
 439     double ix = points[i - 1].x;
 440     double iy = points[i - 1].y;
 441     double iz = points[i - 1].z;
 442     double ox = points[i].x;
 443     double oy = points[i].y;
 444     double oz = points[i].z;
 445     double t[6];
 446     double swap;
 447     double x, y, z;             /* possible intersection point */
 448     double t_min, t_max;
 449
 450     /* nasty degenerate cases, effectively drawing to an infinity point (?)
 451        cope with them here, so don't process them as a "real" OUTRANGE point
 452
 453        If more than one coord is -VERYLARGE, then can't ratio the "infinities"
 454        so drop out by returning FALSE */
 455
 456     count = 0;
 457     if (ix == -VERYLARGE)
 458         count++;
 459     if (ox == -VERYLARGE)
 460         count++;
 461     if (iy == -VERYLARGE)
 462         count++;
 463     if (oy == -VERYLARGE)
 464         count++;
 465     if (iz == -VERYLARGE)
 466         count++;
 467     if (oz == -VERYLARGE)
 468         count++;
 469
 470     /* either doesn't pass through 3D volume *or*
 471        can't ratio infinities to get a direction to draw line, so simply return(FALSE) */
 472     if (count > 1) {
 473         return (FALSE);
 474     }
 475
 476     if (ox == -VERYLARGE || ix == -VERYLARGE) {
 477         if (ix == -VERYLARGE) {
 478             /* swap points so ix/iy/iz don't have a -VERYLARGE component */
 479             x = ix;
 480             ix = ox;
 481             ox = x;
 482             y = iy;
 483             iy = oy;
 484             oy = y;
 485             z = iz;
 486             iz = oz;
 487             oz = z;
 488         }
 489         /* check actually passes through the 3D graph volume */
 490
 491         if (ix > axis_array[FIRST_X_AXIS].max
 492             && IN_AXIS_RANGE(iy, FIRST_Y_AXIS)
 493             && IN_AXIS_RANGE(iz, FIRST_Z_AXIS)) {
 494             lx[0] = axis_array[FIRST_X_AXIS].min;
 495             ly[0] = iy;
 496             lz[0] = iz;
 497
 498             lx[1] = axis_array[FIRST_X_AXIS].max;
 499             ly[1] = iy;
 500             lz[1] = iz;
 501
 502             return (TRUE);
 503         } else {
 504             return (FALSE);
 505         }
 506     }
 507     if (oy == -VERYLARGE || iy == -VERYLARGE) {
 508         if (iy == -VERYLARGE) {
 509             /* swap points so ix/iy/iz don't have a -VERYLARGE component */
 510             x = ix;
 511             ix = ox;
 512             ox = x;
 513             y = iy;
 514             iy = oy;
 515             oy = y;
 516             z = iz;
 517             iz = oz;
 518             oz = z;
 519         }
 520         /* check actually passes through the 3D graph volume */
 521         if (iy > axis_array[FIRST_Y_AXIS].max
 522             && IN_AXIS_RANGE(ix, FIRST_X_AXIS)
 523             && IN_AXIS_RANGE(iz, FIRST_Z_AXIS)) {
 524             lx[0] = ix;
 525             ly[0] = axis_array[FIRST_Y_AXIS].min;
 526             lz[0] = iz;
 527
 528             lx[1] = ix;
 529             ly[1] = axis_array[FIRST_Y_AXIS].max;
 530             lz[1] = iz;
 531
 532             return (TRUE);
 533         } else {
 534             return (FALSE);
 535         }
 536     }
 537     if (oz == -VERYLARGE || iz == -VERYLARGE) {
 538         if (iz == -VERYLARGE) {
 539             /* swap points so ix/iy/iz don't have a -VERYLARGE component */
 540             x = ix;
 541             ix = ox;
 542             ox = x;
 543             y = iy;
 544             iy = oy;
 545             oy = y;
 546             z = iz;
 547             iz = oz;
 548             oz = z;
 549         }
 550         /* check actually passes through the 3D graph volume */
 551         if (iz > axis_array[FIRST_Z_AXIS].max
 552             && IN_AXIS_RANGE(ix, FIRST_X_AXIS)
 553             && IN_AXIS_RANGE(iy, FIRST_Y_AXIS)) {
 554             lx[0] = ix;
 555             ly[0] = iy;
 556             lz[0] = axis_array[FIRST_Z_AXIS].min;
 557
 558             lx[1] = ix;
 559             ly[1] = iy;
 560             lz[1] = axis_array[FIRST_Z_AXIS].max;
 561
 562             return (TRUE);
 563         } else {
 564             return (FALSE);
 565         }
 566     }
 567     /*
 568      * Quick outcode tests on the 3d graph volume
 569      */
 570
 571     /* test z coord first --- most surface OUTRANGE points generated
 572      * between axis_array[FIRST_Z_AXIS].min and baseplane (i.e. when
 573      * ticslevel is non-zero)
 574      */
 575     if (GPMAX(iz, oz) < axis_array[FIRST_Z_AXIS].min
 576         || GPMIN(iz, oz) > axis_array[FIRST_Z_AXIS].max)
 577         return (FALSE);
 578
 579     if (GPMAX(ix, ox) < axis_array[FIRST_X_AXIS].min
 580         || GPMIN(ix, ox) > axis_array[FIRST_X_AXIS].max)
 581         return (FALSE);
 582
 583     if (GPMAX(iy, oy) < axis_array[FIRST_Y_AXIS].min
 584         || GPMIN(iy, oy) > axis_array[FIRST_Y_AXIS].max)
 585         return (FALSE);
 586
 587     /* Special horizontal/vertical, etc. cases are checked and
 588      * remaining slant lines are checked separately.
 589      *
 590      * The slant line intersections are solved using the parametric
 591      * form of the equation for a line, since if we test x/y/z min/max
 592      * planes explicitly then e.g. a line passing through a corner
 593      * point (x_min,y_min,z_min) actually intersects all 3 planes and
 594      * hence further tests would be required to anticipate this and
 595      * similar situations. */
 596
 597     /* Can have case (ix == ox && iy == oy && iz == oz) as both points
 598      * OUTRANGE */
 599     if (ix == ox && iy == oy && iz == oz) {
 600         /* but as only define single outrange point, can't intersect
 601          * 3D graph volume */
 602         return (FALSE);
 603     }
 604
 605     if (ix == ox) {
 606         if (iy == oy) {
 607             /* line parallel to z axis */
 608
 609             /* x and y coords must be in range, and line must span
 610              * both FIRST_Z_AXIS->min and ->max.
 611              *
 612              * note that spanning FIRST_Z_AXIS->min implies spanning
 613              * ->max as both points OUTRANGE */
 614
 615             if (!IN_AXIS_RANGE(ix, FIRST_X_AXIS)
 616                 || !IN_AXIS_RANGE(iy, FIRST_Y_AXIS)) {
 617                 return (FALSE);
 618             }
 619             if (inrange(axis_array[FIRST_Z_AXIS].min, iz, oz)) {
 620                 lx[0] = ix;
 621                 ly[0] = iy;
 622                 lz[0] = axis_array[FIRST_Z_AXIS].min;
 623
 624                 lx[1] = ix;
 625                 ly[1] = iy;
 626                 lz[1] = axis_array[FIRST_Z_AXIS].max;
 627
 628                 return (TRUE);
 629             } else
 630                 return (FALSE);
 631         }
 632         if (iz == oz) {
 633             /* line parallel to y axis */
 634             if (!IN_AXIS_RANGE(ix, FIRST_X_AXIS)
 635                 || !IN_AXIS_RANGE(iz, FIRST_Z_AXIS)) {
 636                 return (FALSE);
 637             }
 638             if (inrange(axis_array[FIRST_Y_AXIS].min, iy, oy)) {
 639                 lx[0] = ix;
 640                 ly[0] = axis_array[FIRST_Y_AXIS].min;
 641                 lz[0] = iz;
 642
 643                 lx[1] = ix;
 644                 ly[1] = axis_array[FIRST_Y_AXIS].max;
 645                 lz[1] = iz;
 646
 647                 return (TRUE);
 648             } else
 649                 return (FALSE);
 650         }
 651
 652
 653         /* nasty 2D slanted line in a yz plane */
 654         if (!IN_AXIS_RANGE(ox, FIRST_X_AXIS))
 655             return (FALSE);
 656
 657         t[0] = (axis_array[FIRST_Y_AXIS].min - iy) / (oy - iy);
 658         t[1] = (axis_array[FIRST_Y_AXIS].max - iy) / (oy - iy);
 659
 660         if (t[0] > t[1]) {
 661             swap = t[0];
 662             t[0] = t[1];
 663             t[1] = swap;
 664         }
 665         t[2] = (axis_array[FIRST_Z_AXIS].min - iz) / (oz - iz);
 666         t[3] = (axis_array[FIRST_Z_AXIS].max - iz) / (oz - iz);
 667
 668         if (t[2] > t[3]) {
 669             swap = t[2];
 670             t[2] = t[3];
 671             t[3] = swap;
 672         }
 673         t_min = GPMAX(GPMAX(t[0], t[2]), 0.0);
 674         t_max = GPMIN(GPMIN(t[1], t[3]), 1.0);
 675
 676         if (t_min > t_max)
 677             return (FALSE);
 678
 679         lx[0] = ix;
 680         ly[0] = iy + t_min * (oy - iy);
 681         lz[0] = iz + t_min * (oz - iz);
 682
 683         lx[1] = ix;
 684         ly[1] = iy + t_max * (oy - iy);
 685         lz[1] = iz + t_max * (oz - iz);
 686
 687         /* Can only have 0 or 2 intersection points -- only need test
 688          * one coord */
 689         if (IN_AXIS_RANGE(ly[0], FIRST_Y_AXIS)
 690             && IN_AXIS_RANGE(lz[0], FIRST_Z_AXIS)) {
 691             return (TRUE);
 692         }
 693         return (FALSE);
 694     }
 695
 696     if (iy == oy) {
 697         /* already checked case (ix == ox && iy == oy) */
 698         if (oz == iz) {
 699             /* line parallel to x axis */
 700             if (!IN_AXIS_RANGE(iy, FIRST_Y_AXIS)
 701                 || !IN_AXIS_RANGE(iz, FIRST_Z_AXIS)) {
 702                 return (FALSE);
 703             }
 704             if (inrange(axis_array[FIRST_X_AXIS].min, ix, ox)) {
 705                 lx[0] = axis_array[FIRST_X_AXIS].min;
 706                 ly[0] = iy;
 707                 lz[0] = iz;
 708
 709                 lx[1] = axis_array[FIRST_X_AXIS].max;
 710                 ly[1] = iy;
 711                 lz[1] = iz;
 712
 713                 return (TRUE);
 714             } else
 715                 return (FALSE);
 716         }
 717         /* nasty 2D slanted line in an xz plane */
 718
 719         if (!IN_AXIS_RANGE(oy, FIRST_Y_AXIS))
 720             return (FALSE);
 721
 722         t[0] = (axis_array[FIRST_X_AXIS].min - ix) / (ox - ix);
 723         t[1] = (axis_array[FIRST_X_AXIS].max - ix) / (ox - ix);
 724
 725         if (t[0] > t[1]) {
 726             swap = t[0];
 727             t[0] = t[1];
 728             t[1] = swap;
 729         }
 730         t[2] = (axis_array[FIRST_Z_AXIS].min - iz) / (oz - iz);
 731         t[3] = (axis_array[FIRST_Z_AXIS].max - iz) / (oz - iz);
 732
 733         if (t[2] > t[3]) {
 734             swap = t[2];
 735             t[2] = t[3];
 736             t[3] = swap;
 737         }
 738         t_min = GPMAX(GPMAX(t[0], t[2]), 0.0);
 739         t_max = GPMIN(GPMIN(t[1], t[3]), 1.0);
 740
 741         if (t_min > t_max)
 742             return (FALSE);
 743
 744         lx[0] = ix + t_min * (ox - ix);
 745         ly[0] = iy;
 746         lz[0] = iz + t_min * (oz - iz);
 747
 748         lx[1] = ix + t_max * (ox - ix);
 749         ly[1] = iy;
 750         lz[1] = iz + t_max * (oz - iz);
 751
 752         /*
 753          * Can only have 0 or 2 intersection points -- only need test one coord
 754          */
 755         if (IN_AXIS_RANGE(lx[0], FIRST_X_AXIS)
 756             && IN_AXIS_RANGE(lz[0], FIRST_Z_AXIS)) {
 757             return (TRUE);
 758         }
 759         return (FALSE);
 760     }
 761     if (iz == oz) {
 762         /* already checked cases (ix == ox && iz == oz) and (iy == oy
 763            && iz == oz) */
 764
 765         /* nasty 2D slanted line in an xy plane */
 766
 767         if (!IN_AXIS_RANGE(oz, FIRST_Z_AXIS))
 768             return (FALSE);
 769
 770         t[0] = (axis_array[FIRST_X_AXIS].min - ix) / (ox - ix);
 771         t[1] = (axis_array[FIRST_X_AXIS].max - ix) / (ox - ix);
 772
 773         if (t[0] > t[1]) {
 774             swap = t[0];
 775             t[0] = t[1];
 776             t[1] = swap;
 777         }
 778         t[2] = (axis_array[FIRST_Y_AXIS].min - iy) / (oy - iy);
 779         t[3] = (axis_array[FIRST_Y_AXIS].max - iy) / (oy - iy);
 780
 781         if (t[2] > t[3]) {
 782             swap = t[2];
 783             t[2] = t[3];
 784             t[3] = swap;
 785         }
 786         t_min = GPMAX(GPMAX(t[0], t[2]), 0.0);
 787         t_max = GPMIN(GPMIN(t[1], t[3]), 1.0);
 788
 789         if (t_min > t_max)
 790             return (FALSE);
 791
 792         lx[0] = ix + t_min * (ox - ix);
 793         ly[0] = iy + t_min * (oy - iy);
 794         lz[0] = iz;
 795
 796         lx[1] = ix + t_max * (ox - ix);
 797         ly[1] = iy + t_max * (oy - iy);
 798         lz[1] = iz;
 799
 800         /*
 801          * Can only have 0 or 2 intersection points -- only need test one coord
 802          */
 803         if (IN_AXIS_RANGE(lx[0], FIRST_X_AXIS)
 804             && IN_AXIS_RANGE(ly[0], FIRST_Y_AXIS)) {
 805             return (TRUE);
 806         }
 807         return (FALSE);
 808     }
 809     /* really nasty general slanted 3D case */
 810
 811     /*
 812        Solve parametric equation
 813
 814        (ix, iy, iz) + t (diff_x, diff_y, diff_z)
 815
 816        where 0.0 <= t <= 1.0 and
 817
 818        diff_x = (ox - ix);
 819        diff_y = (oy - iy);
 820        diff_z = (oz - iz);
 821      */
 822
 823     t[0] = (axis_array[FIRST_X_AXIS].min - ix) / (ox - ix);
 824     t[1] = (axis_array[FIRST_X_AXIS].max - ix) / (ox - ix);
 825
 826     if (t[0] > t[1]) {
 827         swap = t[0];
 828         t[0] = t[1];
 829         t[1] = swap;
 830     }
 831     t[2] = (axis_array[FIRST_Y_AXIS].min - iy) / (oy - iy);
 832     t[3] = (axis_array[FIRST_Y_AXIS].max - iy) / (oy - iy);
 833
 834     if (t[2] > t[3]) {
 835         swap = t[2];
 836         t[2] = t[3];
 837         t[3] = swap;
 838     }
 839     t[4] = (iz == oz) ? 0.0 : (axis_array[FIRST_Z_AXIS].min - iz) / (oz - iz);
 840     t[5] = (iz == oz) ? 1.0 : (axis_array[FIRST_Z_AXIS].max - iz) / (oz - iz);
 841
 842     if (t[4] > t[5]) {
 843         swap = t[4];
 844         t[4] = t[5];
 845         t[5] = swap;
 846     }
 847     t_min = GPMAX(GPMAX(t[0], t[2]), GPMAX(t[4], 0.0));
 848     t_max = GPMIN(GPMIN(t[1], t[3]), GPMIN(t[5], 1.0));
 849
 850     if (t_min > t_max)
 851         return (FALSE);
 852
 853     lx[0] = ix + t_min * (ox - ix);
 854     ly[0] = iy + t_min * (oy - iy);
 855     lz[0] = iz + t_min * (oz - iz);
 856
 857     lx[1] = ix + t_max * (ox - ix);
 858     ly[1] = iy + t_max * (oy - iy);
 859     lz[1] = iz + t_max * (oz - iz);
 860
 861     /*
 862      * Can only have 0 or 2 intersection points -- only need test one coord
 863      */
 864     if (IN_AXIS_RANGE(lx[0], FIRST_X_AXIS)
 865         && IN_AXIS_RANGE(ly[0], FIRST_Y_AXIS)
 866         && IN_AXIS_RANGE(lz[0], FIRST_Z_AXIS)) {
 867         return (TRUE);
 868     }
 869     return (FALSE);
 870 }
 871
 872 /* Performs transformation from 'user coordinates' to a normalized
 873  * vector in 'graph coordinates' (-1..1 in all three directions).  */
 874 void
 875 map3d_xyz(
 876     double x, double y, double z,               /* user coordinates */
 877     p_vertex out)
 878 {
 879     int i, j;
 880     double V[4], Res[4];        /* Homogeneous coords. vectors. */
 881
 882     /* Normalize object space to -1..1 */
 883     V[0] = map_x3d(x);
 884     V[1] = map_y3d(y);
 885     V[2] = map_z3d(z);
 886     V[3] = 1.0;
 887
 888     /* Res[] = V[] * trans_mat[][] (uses row-vectors) */
 889     for (i = 0; i < 4; i++) {
 890         Res[i] = trans_mat[3][i];               /* V[3] is 1. anyway */
 891         for (j = 0; j < 3; j++)
 892             Res[i] += V[j] * trans_mat[j][i];
 893     }
 894
 895     if (Res[3] == 0)
 896         Res[3] = 1.0e-5;
 897
 898     out->x = Res[0] / Res[3];
 899     out->y = Res[1] / Res[3];
 900     out->z = Res[2] / Res[3];
 901     /* store z for later color calculation */
 902     out->real_z = z;
 903 #ifdef EAM_DATASTRINGS
 904     out->label = NULL;
 905 #endif
 906 }
 907
 908
 909 /* DJS (20 Aug 2004):  A more precise double version of map3d_xy() is
 910  * is required for the image routine.  The original intention was to
 911  * reuse the double version of the routine to generate the unsigned
 912  * int versin of the routine.  However, that caused rounding problems
 913  * such that PostScript versions of all the demos didn't come out
 914  * quite exactly the same.
 915  *
 916  * The define switch below will allow either code reuse or code
 917  * replication.  My advice is to study the rounding problem and
 918  * decide if the code reuse rounding is just as well as the code
 919  * replication approach.  If so, go the code reuse route and toss
 920  * the replicated code.
 921  */
 922 #define REPLICATE_CODE_FOR_BACKWARD_COMPATIBLE_ROUNDING 1
 923
 924 #if REPLICATE_CODE_FOR_BACKWARD_COMPATIBLE_ROUNDING
 925
 926 /* Function to map from user 3D space to normalized 'camera' view
 927  * space, and from there directly to terminal coordinates */
 928 void
 929 map3d_xy(
 930     double x, double y, double z,
 931     unsigned int *xt, unsigned int *yt)
 932 {
 933     int i, j;
 934     double v[4], res[4],        /* Homogeneous coords. vectors. */
 935      w = trans_mat[3][3];
 936
 937     v[0] = map_x3d(x);          /* Normalize object space to -1..1 */
 938     v[1] = map_y3d(y);
 939     v[2] = map_z3d(z);
 940     v[3] = 1.0;
 941
 942     for (i = 0; i < 2; i++) {   /* Dont use the third axes (z). */
 943         res[i] = trans_mat[3][i];       /* Initiate it with the weight factor */
 944         for (j = 0; j < 3; j++)
 945             res[i] += v[j] * trans_mat[j][i];
 946     }
 947
 948     for (i = 0; i < 3; i++)
 949         w += v[i] * trans_mat[i][3];
 950     if (w == 0)
 951         w = 1e-5;
 952
 953     if (lmargin.scalex == screen || rmargin.scalex == screen)
 954         *xt = res[0] * xscaler/w + xmiddle;
 955     else
 956         *xt = (unsigned int) ((res[0] * xscaler / w) + xmiddle);
 957
 958     if (tmargin.scalex == screen || bmargin.scalex == screen)
 959         *yt = res[1] * yscaler/w + ymiddle;
 960     else
 961         *yt = (unsigned int) ((res[1] * yscaler / w) + ymiddle);
 962 }
 963
 964 /* Function to map from user 3D space to normalized 'camera' view
 965  * space, and from there directly to terminal coordinates */
 966 void
 967 map3d_xy_double(
 968     double x, double y, double z,
 969     double *xt, double *yt)
 970 {
 971     int i, j;
 972     double v[4], res[4],        /* Homogeneous coords. vectors. */
 973      w = trans_mat[3][3];
 974
 975     v[0] = map_x3d(x);          /* Normalize object space to -1..1 */
 976     v[1] = map_y3d(y);
 977     v[2] = map_z3d(z);
 978     v[3] = 1.0;
 979
 980     for (i = 0; i < 2; i++) {   /* Dont use the third axes (z). */
 981         res[i] = trans_mat[3][i];       /* Initiate it with the weight factor */
 982         for (j = 0; j < 3; j++)
 983             res[i] += v[j] * trans_mat[j][i];
 984     }
 985
 986     for (i = 0; i < 3; i++)
 987         w += v[i] * trans_mat[i][3];
 988     if (w == 0)
 989         w = 1e-5;
 990
 991     if (lmargin.scalex == screen || rmargin.scalex == screen)
 992         *xt = res[0] * xscaler + xmiddle;
 993     else
 994         *xt = (res[0] * xscaler / w) + xmiddle;
 995
 996     if (tmargin.scalex == screen || bmargin.scalex == screen)
 997         *yt = res[1] * yscaler + ymiddle;
 998     else
 999         *yt = (res[1] * yscaler / w) + ymiddle;
1000 }
1001
1002 #else /* REPLICATE_CODE_FOR_BACKWARD_COMPATIBLE_ROUNDING */
1003
1004 /* Function to map from user 3D space to normalized 'camera' view
1005  * space, and from there directly to terminal coordinates */
1006 void
1007 map3d_xy(
1008     double x, double y, double z,
1009     unsigned int *xt, unsigned int *yt)
1010 {
1011 #ifdef WITH_IMAGE
1012     double xtd, ytd;
1013     map3d_xy_double(x, y, z, &xtd, &ytd);
1014     *xt = (unsigned int) xtd;
1015     *yt = (unsigned int) ytd;
1016 }
1017
1018 void
1019 map3d_xy_double(
1020     double x, double y, double z,
1021     double *xt, double *yt)
1022 {
1023 #endif
1024     int i, j;
1025     double v[4], res[4],        /* Homogeneous coords. vectors. */
1026      w = trans_mat[3][3];
1027
1028     v[0] = map_x3d(x);          /* Normalize object space to -1..1 */
1029     v[1] = map_y3d(y);
1030     v[2] = map_z3d(z);
1031     v[3] = 1.0;
1032
1033     for (i = 0; i < 2; i++) {   /* Dont use the third axes (z). */
1034         res[i] = trans_mat[3][i];       /* Initiate it with the weight factor */
1035         for (j = 0; j < 3; j++)
1036             res[i] += v[j] * trans_mat[j][i];
1037     }
1038
1039     for (i = 0; i < 3; i++)
1040         w += v[i] * trans_mat[i][3];
1041     if (w == 0)
1042         w = 1e-5;
1043
1044 #ifdef WITH_IMAGE
1045     *xt = ((res[0] * xscaler / w) + xmiddle);
1046     *yt = ((res[1] * yscaler / w) + ymiddle);
1047 #else
1048     *xt = (unsigned int) ((res[0] * xscaler / w) + xmiddle);
1049     *yt = (unsigned int) ((res[1] * yscaler / w) + ymiddle);
1050 #endif
1051 }
1052
1053 #endif /* REPLICATE_CODE_FOR_BACKWARD_COMPATIBLE_ROUNDING */
1054
1055
1056 /* HBB 20020313: New routine, broken out of draw3d_point, to be used
1057  * to output a single point without any checks for hidden3d */
1058 static GP_INLINE void
1059 draw3d_point_unconditional(p_vertex v, struct lp_style_type *lp)
1060 {
1061     unsigned int x, y;
1062
1063     TERMCOORD(v, x, y);
1064     term_apply_lp_properties(lp);
1065     /* HBB 20010822: implemented "linetype palette" for points, too */
1066     if (lp->use_palette) {
1067         set_color(cb2gray( z2cb(v->real_z) ));
1068     }
1069     if (!clip_point(x, y))
1070         (term->point) (x, y, lp->p_type);
1071 }
1072
1073 /* Moved this upward, to make optional inlining in draw3d_line easier
1074  * for compilers */
1075 /* HBB 20021128: removed GP_INLINE qualifier to avoid MSVC++ silliness */
1076 void
1077 draw3d_line_unconditional(
1078     p_vertex v1, p_vertex v2,
1079     struct lp_style_type *lp,
1080     int linetype)
1081 {
1082     unsigned int x1, y1, x2, y2;
1083     struct lp_style_type ls = *lp;
1084
1085     /* HBB 20020312: v2 can be NULL, if this call is coming from
1086     draw_line_hidden. --> redirect to point drawing routine */
1087     if (! v2) {
1088         draw3d_point_unconditional(v1, lp);
1089         return;
1090     }
1091
1092     TERMCOORD(v1, x1, y1);
1093     TERMCOORD(v2, x2, y2);
1094
1095     /* User-specified line styles */
1096     if (prefer_line_styles && linetype >= 0)
1097         lp_use_properties(&ls, linetype+1, FALSE);
1098
1099     /* The usual case of auto-generated line types */
1100     else
1101         ls.l_type = linetype;
1102
1103     /* Color by Z value */
1104     if (ls.pm3d_color.type == TC_Z)
1105             ls.pm3d_color.value = (v1->real_z + v2->real_z) * 0.5;
1106
1107     term_apply_lp_properties(&ls);
1108     draw_clip_line(x1,y1,x2,y2);
1109 }
1110
1111 void
1112 draw3d_line (p_vertex v1, p_vertex v2, struct lp_style_type *lp)
1113 {
1114 #ifndef LITE
1115     /* hidden3d routine can't work if no surface was drawn at all */
1116     if (hidden3d && draw_surface) {
1117         draw_line_hidden(v1, v2, lp);
1118         return;
1119     }
1120 #endif
1121
1122     draw3d_line_unconditional(v1, v2, lp, lp->l_type);
1123
1124 }
1125
1126 /* HBB 20000621: new routine, to allow for hiding point symbols behind
1127  * the surface */
1128 void
1129 draw3d_point(p_vertex v, struct lp_style_type *lp)
1130 {
1131 #ifndef LITE
1132     /* hidden3d routine can't work if no surface was drawn at all */
1133     if (hidden3d && draw_surface) {
1134         /* Draw vertex as a zero-length edge */
1135         draw_line_hidden(v, NULL, lp);
1136         return;
1137     }
1138 #endif
1139
1140     draw3d_point_unconditional(v, lp);
1141 }
1142
1143 /* HBB NEW 20031218: tools for drawing polylines in 3D with a semantic
1144  * like term->move() and term->vector() */
1145
1146 /* Previous points 3D position */
1147 static vertex polyline3d_previous_vertex;
1148
1149 void
1150 polyline3d_start(p_vertex v1)
1151 {
1152     unsigned int x1, y1;
1153
1154     polyline3d_previous_vertex = *v1;
1155 #ifndef LITE
1156     if (hidden3d && draw_surface)
1157         return;
1158 #endif /* LITE */
1159
1160     TERMCOORD(v1, x1, y1);
1161     /* HBB FIXME 20031219: no clipping!? */
1162     term->move(x1, y1);
1163 }
1164
1165 void
1166 polyline3d_next(p_vertex v2, struct lp_style_type *lp)
1167 {
1168     unsigned int x2, y2;
1169
1170     /* Copied from draw3d_line(): */
1171 #ifndef LITE
1172     /* FIXME HBB 20031218: hidden3d mode will still create isolated
1173      * edges! */
1174     if (hidden3d && draw_surface) {
1175         draw_line_hidden(&polyline3d_previous_vertex, v2, lp);
1176         polyline3d_previous_vertex = *v2;
1177         return;
1178     }
1179 #endif
1180
1181     /* Copied from draw3d_line_unconditional: */
1182     /* If use_palette is active, polylines can't be used -->
1183      * revert back to old method */
1184     if (lp->use_palette) {
1185         draw3d_line_unconditional(&polyline3d_previous_vertex, v2,
1186                                   lp, lp->l_type);
1187         polyline3d_previous_vertex = *v2;
1188         return;
1189
1190     }
1191
1192     TERMCOORD(v2, x2, y2);
1193     /* FIXME HBB 20031219: no clipping?! */
1194     term->vector(x2, y2);
1195
1196     polyline3d_previous_vertex = *v2;
1197 }