2 static char *RCSid() { return RCSid("$Id: util3d.c,v 1.27.2.4 2009/01/07 22:58:27 sfeam Exp $"); }
5 /* GNUPLOT - util3d.c */
8 * Copyright 1986 - 1993, 1998, 2004 Thomas Williams, Colin Kelley
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.
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,
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
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.
33 * This software is provided "as is" without express or implied warranty
34 * to the extent permitted by applicable law.
39 * 19 September 1992 Lawrence Crowl (crowl@cs.orst.edu)
40 * Added user-specified bases for log scaling.
42 * 3.6 - split graph3d.c into graph3d.c (graph),
43 * util3d.c (intersections, etc)
44 * hidden3d.c (hidden-line removal code)
55 /* HBB 990826: all that stuff referenced from other modules is now
56 * exported in graph3d.h, instead of being listed here */
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 *));
63 mat_unit(transform_matrix mat)
67 for (i = 0; i < 4; i++)
68 for (j = 0; j < 4; j++)
75 #if 0 /* HBB 990829: unused --> commented out */
77 mat_trans(double tx, double ty, double tz, transform_matrix mat)
79 mat_unit(mat); /* Make it unit matrix. */
84 #endif /* commented out */
87 mat_scale(double sx, double sy, double sz, transform_matrix mat)
89 mat_unit(mat); /* Make it unit matrix. */
96 mat_rot_x(double teta, transform_matrix mat)
98 double cos_teta, sin_teta;
101 cos_teta = cos(teta);
102 sin_teta = sin(teta);
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;
111 #if 0 /* HBB 990829: unused --> commented out */
113 mat_rot_y(double teta, transform_matrix mat)
115 double cos_teta, sin_teta;
118 cos_teta = cos(teta);
119 sin_teta = sin(teta);
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;
127 #endif /* commented out */
130 mat_rot_z(double teta, transform_matrix mat)
132 double cos_teta, sin_teta;
135 cos_teta = cos(teta);
136 sin_teta = sin(teta);
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;
145 /* Multiply two transform_matrix. Result can be one of two operands. */
148 transform_matrix mat_res,
149 transform_matrix mat1, transform_matrix mat2)
152 transform_matrix mat_res_temp;
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];
160 for (i = 0; i < 4; i++)
161 for (j = 0; j < 4; j++)
162 mat_res[i][j] = mat_res_temp[i][j];
165 #define IN_AXIS_RANGE(val, axis) \
166 inrange((val), axis_array[axis].min, axis_array[axis].max)
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
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 */
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 */
189 if (points[i].type == INRANGE) {
190 /* swap points around so that ix/ix/iz are INRANGE and ox/oy/oz are OUTRANGE */
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
205 If more than one coord is -VERYLARGE, then can't ratio the "infinities"
206 so drop out by returning FALSE */
209 if (ox == -VERYLARGE)
211 if (oy == -VERYLARGE)
213 if (oz == -VERYLARGE)
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 */
230 if (ox == -VERYLARGE) {
231 *ex = AXIS_ACTUAL_MIN(FIRST_X_AXIS);
234 if (oy == -VERYLARGE) {
235 *ey = AXIS_ACTUAL_MIN(FIRST_Y_AXIS);
238 /* obviously oz is -VERYLARGE and (ox != -VERYLARGE && oy != -VERYLARGE) */
239 *ez = AXIS_ACTUAL_MIN(FIRST_Z_AXIS);
243 * Can't have case (ix == ox && iy == oy && iz == oz) as one point
244 * is INRANGE and one point is OUTRANGE.
248 /* line parallel to z axis */
250 /* assume iy in yrange, && ix in xrange */
251 *ex = ix; /* == ox */
252 *ey = iy; /* == oy */
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);
259 graph_error("error in edge3d_intersect");
265 /* line parallel to y axis */
267 /* assume iz in zrange && ix in xrange */
268 *ex = ix; /* == ox */
269 *ez = iz; /* == oz */
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);
276 graph_error("error in edge3d_intersect");
282 /* nasty 2D slanted line in a yz plane */
285 #define INTERSECT_PLANE(cut, axis, eff, eff_axis, res_x, res_y, res_z) \
287 if (inrange(cut, i##axis, o##axis) \
289 && cut != o##axis) { \
290 eff = (cut - i##axis) \
291 * ((o##eff - i##eff) / (o##axis - i##axis)) \
293 if (IN_AXIS_RANGE(eff, eff_axis)) { \
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) */
313 /* already checked case (ix == ox && iy == oy) */
315 /* line parallel to x axis */
317 /* assume inrange(iz) && inrange(iy) */
318 *ey = iy; /* == oy */
319 *ez = iz; /* == oz */
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);
326 graph_error("error in edge3d_intersect");
331 /* nasty 2D slanted line in an xz plane */
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));
344 /* already checked cases (ix == ox && iz == oz) and (iy == oy
347 /* 2D slanted line in an xy plane */
349 /* assume inrange(oz) */
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);
360 #undef INTERSECT_PLANE
362 /* really nasty general slanted 3D case */
364 #define INTERSECT_DIAG(cut, axis, eff, eff_axis, eff2, eff2_axis, \
365 res_x, res_y, res_z) \
367 if (inrange(cut, i##axis, o##axis) \
369 && cut != o##axis) { \
370 eff = (cut - i##axis) \
371 * ((o##eff - i##eff) / (o##axis - i##axis)) \
373 eff2 = (cut - i##axis) \
374 * ((o##eff2 - i##eff2) / (o##axis - i##axis)) \
376 if (IN_AXIS_RANGE(eff, eff_axis) \
377 && IN_AXIS_RANGE(eff2, eff2_axis)) { \
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);
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);
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));
407 #undef INTERSECT_DIAG
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
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).
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 */
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;
447 double x, y, z; /* possible intersection point */
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
453 If more than one coord is -VERYLARGE, then can't ratio the "infinities"
454 so drop out by returning FALSE */
457 if (ix == -VERYLARGE)
459 if (ox == -VERYLARGE)
461 if (iy == -VERYLARGE)
463 if (oy == -VERYLARGE)
465 if (iz == -VERYLARGE)
467 if (oz == -VERYLARGE)
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) */
476 if (ox == -VERYLARGE || ix == -VERYLARGE) {
477 if (ix == -VERYLARGE) {
478 /* swap points so ix/iy/iz don't have a -VERYLARGE component */
489 /* check actually passes through the 3D graph volume */
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;
498 lx[1] = axis_array[FIRST_X_AXIS].max;
507 if (oy == -VERYLARGE || iy == -VERYLARGE) {
508 if (iy == -VERYLARGE) {
509 /* swap points so ix/iy/iz don't have a -VERYLARGE component */
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)) {
525 ly[0] = axis_array[FIRST_Y_AXIS].min;
529 ly[1] = axis_array[FIRST_Y_AXIS].max;
537 if (oz == -VERYLARGE || iz == -VERYLARGE) {
538 if (iz == -VERYLARGE) {
539 /* swap points so ix/iy/iz don't have a -VERYLARGE component */
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)) {
556 lz[0] = axis_array[FIRST_Z_AXIS].min;
560 lz[1] = axis_array[FIRST_Z_AXIS].max;
568 * Quick outcode tests on the 3d graph volume
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)
575 if (GPMAX(iz, oz) < axis_array[FIRST_Z_AXIS].min
576 || GPMIN(iz, oz) > axis_array[FIRST_Z_AXIS].max)
579 if (GPMAX(ix, ox) < axis_array[FIRST_X_AXIS].min
580 || GPMIN(ix, ox) > axis_array[FIRST_X_AXIS].max)
583 if (GPMAX(iy, oy) < axis_array[FIRST_Y_AXIS].min
584 || GPMIN(iy, oy) > axis_array[FIRST_Y_AXIS].max)
587 /* Special horizontal/vertical, etc. cases are checked and
588 * remaining slant lines are checked separately.
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. */
597 /* Can have case (ix == ox && iy == oy && iz == oz) as both points
599 if (ix == ox && iy == oy && iz == oz) {
600 /* but as only define single outrange point, can't intersect
607 /* line parallel to z axis */
609 /* x and y coords must be in range, and line must span
610 * both FIRST_Z_AXIS->min and ->max.
612 * note that spanning FIRST_Z_AXIS->min implies spanning
613 * ->max as both points OUTRANGE */
615 if (!IN_AXIS_RANGE(ix, FIRST_X_AXIS)
616 || !IN_AXIS_RANGE(iy, FIRST_Y_AXIS)) {
619 if (inrange(axis_array[FIRST_Z_AXIS].min, iz, oz)) {
622 lz[0] = axis_array[FIRST_Z_AXIS].min;
626 lz[1] = axis_array[FIRST_Z_AXIS].max;
633 /* line parallel to y axis */
634 if (!IN_AXIS_RANGE(ix, FIRST_X_AXIS)
635 || !IN_AXIS_RANGE(iz, FIRST_Z_AXIS)) {
638 if (inrange(axis_array[FIRST_Y_AXIS].min, iy, oy)) {
640 ly[0] = axis_array[FIRST_Y_AXIS].min;
644 ly[1] = axis_array[FIRST_Y_AXIS].max;
653 /* nasty 2D slanted line in a yz plane */
654 if (!IN_AXIS_RANGE(ox, FIRST_X_AXIS))
657 t[0] = (axis_array[FIRST_Y_AXIS].min - iy) / (oy - iy);
658 t[1] = (axis_array[FIRST_Y_AXIS].max - iy) / (oy - iy);
665 t[2] = (axis_array[FIRST_Z_AXIS].min - iz) / (oz - iz);
666 t[3] = (axis_array[FIRST_Z_AXIS].max - iz) / (oz - iz);
673 t_min = GPMAX(GPMAX(t[0], t[2]), 0.0);
674 t_max = GPMIN(GPMIN(t[1], t[3]), 1.0);
680 ly[0] = iy + t_min * (oy - iy);
681 lz[0] = iz + t_min * (oz - iz);
684 ly[1] = iy + t_max * (oy - iy);
685 lz[1] = iz + t_max * (oz - iz);
687 /* Can only have 0 or 2 intersection points -- only need test
689 if (IN_AXIS_RANGE(ly[0], FIRST_Y_AXIS)
690 && IN_AXIS_RANGE(lz[0], FIRST_Z_AXIS)) {
697 /* already checked case (ix == ox && iy == oy) */
699 /* line parallel to x axis */
700 if (!IN_AXIS_RANGE(iy, FIRST_Y_AXIS)
701 || !IN_AXIS_RANGE(iz, FIRST_Z_AXIS)) {
704 if (inrange(axis_array[FIRST_X_AXIS].min, ix, ox)) {
705 lx[0] = axis_array[FIRST_X_AXIS].min;
709 lx[1] = axis_array[FIRST_X_AXIS].max;
717 /* nasty 2D slanted line in an xz plane */
719 if (!IN_AXIS_RANGE(oy, FIRST_Y_AXIS))
722 t[0] = (axis_array[FIRST_X_AXIS].min - ix) / (ox - ix);
723 t[1] = (axis_array[FIRST_X_AXIS].max - ix) / (ox - ix);
730 t[2] = (axis_array[FIRST_Z_AXIS].min - iz) / (oz - iz);
731 t[3] = (axis_array[FIRST_Z_AXIS].max - iz) / (oz - iz);
738 t_min = GPMAX(GPMAX(t[0], t[2]), 0.0);
739 t_max = GPMIN(GPMIN(t[1], t[3]), 1.0);
744 lx[0] = ix + t_min * (ox - ix);
746 lz[0] = iz + t_min * (oz - iz);
748 lx[1] = ix + t_max * (ox - ix);
750 lz[1] = iz + t_max * (oz - iz);
753 * Can only have 0 or 2 intersection points -- only need test one coord
755 if (IN_AXIS_RANGE(lx[0], FIRST_X_AXIS)
756 && IN_AXIS_RANGE(lz[0], FIRST_Z_AXIS)) {
762 /* already checked cases (ix == ox && iz == oz) and (iy == oy
765 /* nasty 2D slanted line in an xy plane */
767 if (!IN_AXIS_RANGE(oz, FIRST_Z_AXIS))
770 t[0] = (axis_array[FIRST_X_AXIS].min - ix) / (ox - ix);
771 t[1] = (axis_array[FIRST_X_AXIS].max - ix) / (ox - ix);
778 t[2] = (axis_array[FIRST_Y_AXIS].min - iy) / (oy - iy);
779 t[3] = (axis_array[FIRST_Y_AXIS].max - iy) / (oy - iy);
786 t_min = GPMAX(GPMAX(t[0], t[2]), 0.0);
787 t_max = GPMIN(GPMIN(t[1], t[3]), 1.0);
792 lx[0] = ix + t_min * (ox - ix);
793 ly[0] = iy + t_min * (oy - iy);
796 lx[1] = ix + t_max * (ox - ix);
797 ly[1] = iy + t_max * (oy - iy);
801 * Can only have 0 or 2 intersection points -- only need test one coord
803 if (IN_AXIS_RANGE(lx[0], FIRST_X_AXIS)
804 && IN_AXIS_RANGE(ly[0], FIRST_Y_AXIS)) {
809 /* really nasty general slanted 3D case */
812 Solve parametric equation
814 (ix, iy, iz) + t (diff_x, diff_y, diff_z)
816 where 0.0 <= t <= 1.0 and
823 t[0] = (axis_array[FIRST_X_AXIS].min - ix) / (ox - ix);
824 t[1] = (axis_array[FIRST_X_AXIS].max - ix) / (ox - ix);
831 t[2] = (axis_array[FIRST_Y_AXIS].min - iy) / (oy - iy);
832 t[3] = (axis_array[FIRST_Y_AXIS].max - iy) / (oy - iy);
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);
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));
853 lx[0] = ix + t_min * (ox - ix);
854 ly[0] = iy + t_min * (oy - iy);
855 lz[0] = iz + t_min * (oz - iz);
857 lx[1] = ix + t_max * (ox - ix);
858 ly[1] = iy + t_max * (oy - iy);
859 lz[1] = iz + t_max * (oz - iz);
862 * Can only have 0 or 2 intersection points -- only need test one coord
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)) {
872 /* Performs transformation from 'user coordinates' to a normalized
873 * vector in 'graph coordinates' (-1..1 in all three directions). */
876 double x, double y, double z, /* user coordinates */
880 double V[4], Res[4]; /* Homogeneous coords. vectors. */
882 /* Normalize object space to -1..1 */
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];
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 */
903 #ifdef EAM_DATASTRINGS
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.
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.
922 #define REPLICATE_CODE_FOR_BACKWARD_COMPATIBLE_ROUNDING 1
924 #if REPLICATE_CODE_FOR_BACKWARD_COMPATIBLE_ROUNDING
926 /* Function to map from user 3D space to normalized 'camera' view
927 * space, and from there directly to terminal coordinates */
930 double x, double y, double z,
931 unsigned int *xt, unsigned int *yt)
934 double v[4], res[4], /* Homogeneous coords. vectors. */
937 v[0] = map_x3d(x); /* Normalize object space to -1..1 */
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];
948 for (i = 0; i < 3; i++)
949 w += v[i] * trans_mat[i][3];
953 if (lmargin.scalex == screen || rmargin.scalex == screen)
954 *xt = res[0] * xscaler/w + xmiddle;
956 *xt = (unsigned int) ((res[0] * xscaler / w) + xmiddle);
958 if (tmargin.scalex == screen || bmargin.scalex == screen)
959 *yt = res[1] * yscaler/w + ymiddle;
961 *yt = (unsigned int) ((res[1] * yscaler / w) + ymiddle);
964 /* Function to map from user 3D space to normalized 'camera' view
965 * space, and from there directly to terminal coordinates */
968 double x, double y, double z,
969 double *xt, double *yt)
972 double v[4], res[4], /* Homogeneous coords. vectors. */
975 v[0] = map_x3d(x); /* Normalize object space to -1..1 */
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];
986 for (i = 0; i < 3; i++)
987 w += v[i] * trans_mat[i][3];
991 if (lmargin.scalex == screen || rmargin.scalex == screen)
992 *xt = res[0] * xscaler + xmiddle;
994 *xt = (res[0] * xscaler / w) + xmiddle;
996 if (tmargin.scalex == screen || bmargin.scalex == screen)
997 *yt = res[1] * yscaler + ymiddle;
999 *yt = (res[1] * yscaler / w) + ymiddle;
1002 #else /* REPLICATE_CODE_FOR_BACKWARD_COMPATIBLE_ROUNDING */
1004 /* Function to map from user 3D space to normalized 'camera' view
1005 * space, and from there directly to terminal coordinates */
1008 double x, double y, double z,
1009 unsigned int *xt, unsigned int *yt)
1013 map3d_xy_double(x, y, z, &xtd, &ytd);
1014 *xt = (unsigned int) xtd;
1015 *yt = (unsigned int) ytd;
1020 double x, double y, double z,
1021 double *xt, double *yt)
1025 double v[4], res[4], /* Homogeneous coords. vectors. */
1026 w = trans_mat[3][3];
1028 v[0] = map_x3d(x); /* Normalize object space to -1..1 */
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];
1039 for (i = 0; i < 3; i++)
1040 w += v[i] * trans_mat[i][3];
1045 *xt = ((res[0] * xscaler / w) + xmiddle);
1046 *yt = ((res[1] * yscaler / w) + ymiddle);
1048 *xt = (unsigned int) ((res[0] * xscaler / w) + xmiddle);
1049 *yt = (unsigned int) ((res[1] * yscaler / w) + ymiddle);
1053 #endif /* REPLICATE_CODE_FOR_BACKWARD_COMPATIBLE_ROUNDING */
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)
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) ));
1069 if (!clip_point(x, y))
1070 (term->point) (x, y, lp->p_type);
1073 /* Moved this upward, to make optional inlining in draw3d_line easier
1075 /* HBB 20021128: removed GP_INLINE qualifier to avoid MSVC++ silliness */
1077 draw3d_line_unconditional(
1078 p_vertex v1, p_vertex v2,
1079 struct lp_style_type *lp,
1082 unsigned int x1, y1, x2, y2;
1083 struct lp_style_type ls = *lp;
1085 /* HBB 20020312: v2 can be NULL, if this call is coming from
1086 draw_line_hidden. --> redirect to point drawing routine */
1088 draw3d_point_unconditional(v1, lp);
1092 TERMCOORD(v1, x1, y1);
1093 TERMCOORD(v2, x2, y2);
1095 /* User-specified line styles */
1096 if (prefer_line_styles && linetype >= 0)
1097 lp_use_properties(&ls, linetype+1, FALSE);
1099 /* The usual case of auto-generated line types */
1101 ls.l_type = linetype;
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;
1107 term_apply_lp_properties(&ls);
1108 draw_clip_line(x1,y1,x2,y2);
1112 draw3d_line (p_vertex v1, p_vertex v2, struct lp_style_type *lp)
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);
1122 draw3d_line_unconditional(v1, v2, lp, lp->l_type);
1126 /* HBB 20000621: new routine, to allow for hiding point symbols behind
1129 draw3d_point(p_vertex v, struct lp_style_type *lp)
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);
1140 draw3d_point_unconditional(v, lp);
1143 /* HBB NEW 20031218: tools for drawing polylines in 3D with a semantic
1144 * like term->move() and term->vector() */
1146 /* Previous points 3D position */
1147 static vertex polyline3d_previous_vertex;
1150 polyline3d_start(p_vertex v1)
1152 unsigned int x1, y1;
1154 polyline3d_previous_vertex = *v1;
1156 if (hidden3d && draw_surface)
1160 TERMCOORD(v1, x1, y1);
1161 /* HBB FIXME 20031219: no clipping!? */
1166 polyline3d_next(p_vertex v2, struct lp_style_type *lp)
1168 unsigned int x2, y2;
1170 /* Copied from draw3d_line(): */
1172 /* FIXME HBB 20031218: hidden3d mode will still create isolated
1174 if (hidden3d && draw_surface) {
1175 draw_line_hidden(&polyline3d_previous_vertex, v2, lp);
1176 polyline3d_previous_vertex = *v2;
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,
1187 polyline3d_previous_vertex = *v2;
1192 TERMCOORD(v2, x2, y2);
1193 /* FIXME HBB 20031219: no clipping?! */
1194 term->vector(x2, y2);
1196 polyline3d_previous_vertex = *v2;