2 static char *RCSid() { return RCSid("$Id: pm3d.c,v 1.66.2.6 2009/06/19 00:01:00 sfeam Exp $"); }
9 * Petr Mikulik, since December 1998
10 * Copyright: open source as much as possible
12 * What is here: global variables and routines for the pm3d splotting mode.
13 * This file is included only if PM3D is defined.
25 #include "hidden3d.h" /* p_vertex & map3d_xyz() */
28 #include "setshow.h" /* for surface_rot_z */
29 #include "term_api.h" /* for lp_use_properties() */
30 #include "command.h" /* for c_token */
32 #include <stdlib.h> /* qsort() */
36 Global options for pm3d algorithm (to be accessed by set / show).
41 PM3D_FLUSH_BEGIN, /* flush */
42 0, /* no flushing triangles */
43 PM3D_SCANS_AUTOMATIC, /* scans direction is determined automatically */
44 PM3D_CLIP_4IN, /* clipping: all 4 points in ranges */
45 0, /* no pm3d hidden3d is drawn */
47 0, /* solid (off by default, that means `transparent') */
48 #endif /* PM3D_HAVE_SOLID */
49 PM3D_EXPLICIT, /* implicit */
50 PM3D_WHICHCORNER_MEAN, /* color from which corner(s) */
51 1, /* interpolate along scanline */
52 1 /* interpolate between scanlines */
57 double z; /* maximal z value after rotation to graph coordinate system */
59 gpiPoint icorners[4]; /* also if EXTENDED_COLOR_SPECS is not defined */
62 static int allocated_quadrangles = 0;
63 static int current_quadrangle = 0;
64 static quadrangle* quadrangles = (quadrangle*)0;
66 /* Internal prototypes for this module */
67 static TBOOLEAN plot_has_palette;
68 static double geomean4 __PROTO((double, double, double, double));
69 static double median4 __PROTO((double, double, double, double));
70 static void pm3d_plot __PROTO((struct surface_points *, int));
71 static void pm3d_option_at_error __PROTO((void));
72 static void pm3d_rearrange_part __PROTO((struct iso_curve *, const int, struct iso_curve ***, int *));
73 static void filled_color_contour_plot __PROTO((struct surface_points *, int));
79 /* Geometrical mean = pow( prod(x_i > 0) x_i, 1/N )
80 * Sign of the result: result is positive if 3 or 4 x_i are positive,
81 * it is negative if 3 or all 4 x_i are negative. Helps to splot surface
82 * with all color coordinates negative.
85 geomean4 (double x1, double x2, double x3, double x4)
88 /* return 0 if any of the number is negative */
95 /* honor signess, i.e. sign(geomean) = sign(prod(x_i)) */
96 int neg = (x1 < 0) + (x2 < 0) + (x3 < 0) + (x4 < 0);
98 if (x1 == 0) return 0;
99 /* pow(x, 0.25) is slightly faster than sqrt(sqrt(x)) */
100 x1 = sqrt(sqrt(fabs(x1)));
102 /* such a warning could be helpful, but under normal usage it is just an overhead */
103 if (neg > 1 && interactive && notwarned) {
104 int notwarned = 1; ... to be set on every new splot
106 int_warn(NO_CARET, "corners2color geomean with negative data points");
110 return (neg <= 2) ? x1 : -x1;
115 /* Median: sort values, and then: for N odd, it is the middle value; for N even,
116 * it is mean of the two middle values.
119 median4 (double x1, double x2, double x3, double x4)
122 /* sort them: x1 < x2 and x3 < x4 */
123 if (x1 > x2) { tmp = x2; x2 = x1; x1 = tmp; }
124 if (x3 > x4) { tmp = x3; x3 = x4; x4 = tmp; }
125 /* sum middle numbers */
126 tmp = (x1 < x3) ? x3 : x1;
127 tmp += (x2 < x4) ? x2 : x4;
132 /* Minimum of 4 numbers.
135 minimum4 (double x1, double x2, double x3, double x4)
139 return GPMIN(x1, x3);
143 /* Maximum of 4 numbers.
146 maximum4 (double x1, double x2, double x3, double x4)
150 return GPMAX(x1, x3);
155 * Now the routines which are really just those for pm3d.c
159 * Rescale z to cb values. Nothing to do if both z and cb are linear or log of the
160 * same base, other it has to un-log z and subsequently log it again.
165 if (!Z_AXIS.log && !CB_AXIS.log) /* both are linear */
167 if (Z_AXIS.log && !CB_AXIS.log) /* log z, linear cb */
168 return exp(z * Z_AXIS.log_base); /* unlog(z) */
169 if (!Z_AXIS.log && CB_AXIS.log) /* linear z, log cb */
170 return (log(z) / CB_AXIS.log_base);
172 if (Z_AXIS.base==CB_AXIS.base) /* can we compare double numbers like that? */
174 return z * Z_AXIS.log_base / CB_AXIS.log_base; /* log_cb(unlog_z(z)) */
179 * Rescale cb (color) value into the interval of grays [0,1], taking care
180 * of palette being positive or negative.
181 * Note that it is OK for logarithmic cb-axis too.
186 if (cb <= CB_AXIS.min)
187 return (sm_palette.positive == SMPAL_POSITIVE) ? 0 : 1;
188 if (cb >= CB_AXIS.max)
189 return (sm_palette.positive == SMPAL_POSITIVE) ? 1 : 0;
190 cb = (cb - CB_AXIS.min)
191 / (CB_AXIS.max - CB_AXIS.min);
192 return (sm_palette.positive == SMPAL_POSITIVE) ? cb : 1-cb;
201 struct iso_curve *src,
203 struct iso_curve ***dest,
206 struct iso_curve *scanA;
207 struct iso_curve *scanB;
208 struct iso_curve **scan_array;
210 int invert_order = 0;
212 /* loop over scans in one surface
213 Scans are linked from this_plot->iso_crvs in the opposite order than
214 they are in the datafile.
215 Therefore it is necessary to make vector scan_array of iso_curves.
216 Scans are sorted in scan_array according to pm3d.direction (this can
217 be PM3D_SCANS_FORWARD or PM3D_SCANS_BACKWARD).
219 scan_array = *dest = gp_alloc(len * sizeof(scanA), "pm3d scan array");
221 if (pm3d.direction == PM3D_SCANS_AUTOMATIC) {
224 TBOOLEAN exit_outer_loop = 0;
226 for (scanA = src; scanA && 0 == exit_outer_loop; scanA = scanA->next, len2--) {
231 if ((cnt = scanA->p_count - 1) <= 0)
234 /* ordering within one scan */
235 for (from=0; from<=cnt; from++) /* find 1st non-undefined point */
236 if (scanA->points[from].type != UNDEFINED) {
237 map3d_xyz(scanA->points[from].x, scanA->points[from].y, 0, &vA);
240 for (i=cnt; i>from; i--) /* find the last non-undefined point */
241 if (scanA->points[i].type != UNDEFINED) {
242 map3d_xyz(scanA->points[i].x, scanA->points[i].y, 0, &vA2);
246 if (i - from > cnt * 0.1)
247 /* it is completely arbitrary to request at least
248 * 10% valid samples in this scan. (joze Jun-05-2002) */
249 *invert = (vA2.z > vA.z) ? 0 : 1;
251 continue; /* all points were undefined, so check next scan */
254 /* check the z ordering between scans
255 * Find last scan. If this scan has all points undefined,
256 * find last but one scan, an so on. */
258 for (; len2 >= 3 && !exit_outer_loop; len2--) {
259 for (scanB = scanA-> next, i = len2 - 2; i && scanB; i--)
260 scanB = scanB->next; /* skip over to last scan */
261 if (scanB && scanB->p_count) {
263 for (i = from /* we compare vA.z with vB.z */; i<scanB->p_count; i++) {
264 /* find 1st non-undefined point */
265 if (scanB->points[i].type != UNDEFINED) {
266 map3d_xyz(scanB->points[i].x, scanB->points[i].y, 0, &vB);
267 invert_order = (vB.z > vA.z) ? 0 : 1;
278 fprintf(stderr, "(pm3d_rearrange_part) invert = %d\n", *invert);
279 fprintf(stderr, "(pm3d_rearrange_part) invert_order = %d\n", invert_order);
282 for (scanA = src, scan = len - 1, i = 0; scan >= 0; --scan, i++) {
283 if (pm3d.direction == PM3D_SCANS_AUTOMATIC) {
284 switch (invert_order) {
286 scan_array[scan] = scanA;
290 scan_array[i] = scanA;
293 } else if (pm3d.direction == PM3D_SCANS_FORWARD)
294 scan_array[scan] = scanA;
295 else /* PM3D_SCANS_BACKWARD: i counts scans */
296 scan_array[i] = scanA;
303 * Rearrange scan array
305 * Allocates *first_ptr (and eventually *second_ptr)
306 * which must be freed by the caller
309 pm3d_rearrange_scan_array(
310 struct surface_points *this_plot,
311 struct iso_curve ***first_ptr,
312 int *first_n, int *first_invert,
313 struct iso_curve ***second_ptr,
314 int *second_n, int *second_invert)
317 pm3d_rearrange_part(this_plot->iso_crvs, this_plot->num_iso_read, first_ptr, first_invert);
318 *first_n = this_plot->num_iso_read;
321 struct iso_curve *icrvs = this_plot->iso_crvs;
322 struct iso_curve *icrvs2;
324 /* advance until second part */
325 for (i = 0; i < this_plot->num_iso_read; i++)
327 /* count the number of scans of second part */
328 for (i = 0, icrvs2 = icrvs; icrvs2; icrvs2 = icrvs2->next)
332 pm3d_rearrange_part(icrvs, i, second_ptr, second_invert);
334 *second_ptr = (struct iso_curve **) 0;
339 static int compare_quadrangles(const void* v1, const void* v2)
341 const quadrangle* q1 = (const quadrangle*)v1;
342 const quadrangle* q2 = (const quadrangle*)v2;
346 else if (q1->z < q2->z)
352 void pm3d_depth_queue_clear(void)
354 if (pm3d.direction != PM3D_DEPTH)
359 quadrangles = (quadrangle*)0;
360 allocated_quadrangles = 0;
361 current_quadrangle = 0;
364 void pm3d_depth_queue_flush(void)
366 if (pm3d.direction != PM3D_DEPTH)
369 if (current_quadrangle > 0 && quadrangles) {
376 double z = 0; /* assignment keeps the compiler happy */
377 double w = trans_mat[3][3];
380 for (qp = quadrangles, qe = quadrangles + current_quadrangle; qp != qe; qp++) {
382 gpdPtr = qp->corners;
383 gpiPtr = qp->icorners;
385 for (i = 0; i < 4; i++, gpdPtr++, gpiPtr++) {
387 map3d_xyz(gpdPtr->x, gpdPtr->y, gpdPtr->z, &out);
389 if (i == 0 || out.z > z)
392 gpiPtr->x = (unsigned int) ((out.x * xscaler / w) + xmiddle);
393 gpiPtr->y = (unsigned int) ((out.y * yscaler / w) + ymiddle);
396 qp->z = z; /* maximal z value of all four corners */
399 qsort(quadrangles, current_quadrangle, sizeof (quadrangle), compare_quadrangles);
401 for (qp = quadrangles, qe = quadrangles + current_quadrangle; qp != qe; qp++) {
404 ifilled_quadrangle(qp->icorners);
408 pm3d_depth_queue_clear();
412 * Now the implementation of the pm3d (s)plotting mode
414 * Note: the input parameter at_which_z is char, but an old HP cc requires
415 * ANSI C K&R routines with int only.
418 pm3d_plot(struct surface_points *this_plot, int at_which_z)
420 int j, i, i1, ii, ii1, from, curve, scan, up_to, up_to_minus, invert = 0;
421 int go_over_pts, max_pts;
422 int are_ftriangles, ftriangles_low_pt = -999, ftriangles_high_pt = -999;
423 struct iso_curve *scanA, *scanB;
424 struct coordinate GPHUGE *pointsA, *pointsB;
425 struct iso_curve **scan_array;
428 double cb1, cb2, cb3, cb4;
430 #ifdef EXTENDED_COLOR_SPECS
431 gpiPoint icorners[4];
433 gpdPoint **bl_point = NULL; /* used for bilinear interpolation */
435 /* just a shortcut */
436 TBOOLEAN color_from_column = this_plot->pm3d_color_from_column;
438 if (this_plot == NULL)
441 if (at_which_z != PM3D_AT_BASE && at_which_z != PM3D_AT_TOP && at_which_z != PM3D_AT_SURFACE)
444 /* return if the terminal does not support filled polygons */
445 if (!term->filled_polygon)
448 switch (at_which_z) {
450 corners[0].z = corners[1].z = corners[2].z = corners[3].z = base_z;
453 corners[0].z = corners[1].z = corners[2].z = corners[3].z = ceiling_z;
455 /* the 3rd possibility is surface, PM3D_AT_SURFACE, coded below */
458 scanA = this_plot->iso_crvs;
461 pm3d_rearrange_scan_array(this_plot, &scan_array, &scan_array_n, &invert, (struct iso_curve ***) 0, (int *) 0, (int *) 0);
463 if (pm3d.direction == PM3D_DEPTH) {
465 for (scan = 0; scan < this_plot->num_iso_read - 1; scan++) {
467 scanA = scan_array[scan];
468 scanB = scan_array[scan + 1];
470 are_ftriangles = pm3d.ftriangles && (scanA->p_count != scanB->p_count);
472 allocated_quadrangles += GPMIN(scanA->p_count, scanB->p_count) - 1;
474 allocated_quadrangles += GPMAX(scanA->p_count, scanB->p_count) - 1;
477 allocated_quadrangles *= (pm3d.interp_i > 1) ? pm3d.interp_i : 1;
478 allocated_quadrangles *= (pm3d.interp_j > 1) ? pm3d.interp_j : 1;
479 quadrangles = (quadrangle*)gp_realloc(quadrangles, allocated_quadrangles * sizeof (quadrangle), "pm3d_plot->quadrangles");
480 /* DEBUG: fprintf(stderr, "allocated_quadrangles = %d\n", allocated_quadrangles); */
482 /* pm3d_rearrange_scan_array(this_plot, (struct iso_curve***)0, (int*)0, &scan_array, &invert); */
485 /* debugging: print scan_array */
486 for (scan = 0; scan < this_plot->num_iso_read; scan++) {
487 printf("**** SCAN=%d points=%d\n", scan, scan_array[scan]->p_count);
492 /* debugging: this loop prints properties of all scans */
493 for (scan = 0; scan < this_plot->num_iso_read; scan++) {
494 struct coordinate GPHUGE *points;
495 scanA = scan_array[scan];
496 printf("\n#IsoCurve = scan nb %d, %d points\n#x y z type(in,out,undef)\n", scan, scanA->p_count);
497 for (i = 0, points = scanA->points; i < scanA->p_count; i++) {
498 printf("%g %g %g %c\n",
499 points[i].x, points[i].y, points[i].z, points[i].type == INRANGE ? 'i' : points[i].type == OUTRANGE ? 'o' : 'u');
500 /* Note: INRANGE, OUTRANGE, UNDEFINED */
507 * if bilinear interpolation is enabled, allocate memory for the
508 * interpolated points here
510 if (pm3d.interp_i > 1 || pm3d.interp_j > 1) {
511 bl_point = (gpdPoint **)gp_alloc(sizeof(gpdPoint*) * (pm3d.interp_i+1), "bl-interp along scan");
512 for (i1 = 0; i1 <= pm3d.interp_i; i1++)
513 bl_point[i1] = (gpdPoint *)gp_alloc(sizeof(gpdPoint) * (pm3d.interp_j+1), "bl-interp between scan");
517 * this loop does the pm3d draw of joining two curves
519 * How the loop below works:
520 * - scanB = scan last read; scanA = the previous one
521 * - link the scan from A to B, then move B to A, then read B, then draw
523 for (scan = 0; scan < this_plot->num_iso_read - 1; scan++) {
524 scanA = scan_array[scan];
525 scanB = scan_array[scan + 1];
527 printf("\n#IsoCurveA = scan nb %d has %d points ScanB has %d points\n", scan, scanA->p_count, scanB->p_count);
529 pointsA = scanA->points;
530 pointsB = scanB->points;
531 /* if the number of points in both scans is not the same, then the
532 * starting index (offset) of scan B according to the flushing setting
533 * has to be determined
535 from = 0; /* default is pm3d.flush==PM3D_FLUSH_BEGIN */
536 if (pm3d.flush == PM3D_FLUSH_END)
537 from = abs(scanA->p_count - scanB->p_count);
538 else if (pm3d.flush == PM3D_FLUSH_CENTER)
539 from = abs(scanA->p_count - scanB->p_count) / 2;
540 /* find the minimal number of points in both scans */
541 up_to = GPMIN(scanA->p_count, scanB->p_count) - 1;
542 up_to_minus = up_to - 1; /* calculate only once */
543 are_ftriangles = pm3d.ftriangles && (scanA->p_count != scanB->p_count);
547 max_pts = GPMAX(scanA->p_count, scanB->p_count);
548 go_over_pts = max_pts - 1;
549 /* the j-subrange of quadrangles; in the remaing of the interval
550 * [0..up_to] the flushing triangles are to be drawn */
551 ftriangles_low_pt = from;
552 ftriangles_high_pt = from + up_to_minus;
555 * - the minimal number of points from both scans, if only quadrangles.
556 * - the maximal number of points from both scans if flush triangles
557 * (the missing points in the scan of lower nb of points will be
558 * duplicated from the begin/end points).
560 * Notice: if it would be once necessary to go over points in `backward'
561 * direction, then the loop body below would require to replace the data
562 * point indices `i' by `up_to-i' and `i+1' by `up_to-i-1'.
564 for (j = 0; j < go_over_pts; j++) {
565 /* Now i be the index of the scan with smaller number of points,
566 * ii of the scan with larger number of points. */
567 if (are_ftriangles && (j < ftriangles_low_pt || j > ftriangles_high_pt)) {
568 i = (j <= ftriangles_low_pt) ? 0 : ftriangles_high_pt-from+1;
573 int jj = are_ftriangles ? j - from : j;
575 if (PM3D_SCANS_AUTOMATIC == pm3d.direction && invert)
576 i = up_to_minus - jj;
581 /* From here, i is index to scan A, ii to scan B */
582 if (scanA->p_count > scanB->p_count) {
591 fprintf(stderr,"j=%i: i=%i i1=%i [%i] ii=%i ii1=%i [%i]\n",j,i,i1,scanA->p_count,ii,ii1,scanB->p_count);
593 /* choose the clipping method */
594 if (pm3d.clip == PM3D_CLIP_4IN) {
595 /* (1) all 4 points of the quadrangle must be in x and y range */
596 if (!(pointsA[i].type == INRANGE && pointsA[i1].type == INRANGE &&
597 pointsB[ii].type == INRANGE && pointsB[ii1].type == INRANGE))
599 } else { /* (pm3d.clip == PM3D_CLIP_1IN) */
600 /* (2) all 4 points of the quadrangle must be defined */
601 if (pointsA[i].type == UNDEFINED || pointsA[i1].type == UNDEFINED ||
602 pointsB[ii].type == UNDEFINED || pointsB[ii1].type == UNDEFINED)
604 /* and at least 1 point of the quadrangle must be in x and y range */
605 if (pointsA[i].type == OUTRANGE && pointsA[i1].type == OUTRANGE &&
606 pointsB[ii].type == OUTRANGE && pointsB[ii1].type == OUTRANGE)
610 if ((pm3d.interp_i <= 1 && pm3d.interp_j <= 1) || pm3d.direction == PM3D_DEPTH) {
611 #ifdef EXTENDED_COLOR_SPECS
612 if (!supply_extended_color_specs) {
614 /* Get the gray as the average of the corner z- or gray-positions
615 (note: log scale is already included). The average is calculated here
616 if there is no interpolation (including the "pm3d depthorder" option),
617 otherwise it is done for each interpolated quadrangle later.
618 I always wonder what is faster: d*0.25 or d/4? Someone knows? -- 0.25 (joze) */
619 if (color_from_column) {
620 /* color is set in plot3d.c:get_3ddata() */
621 cb1 = pointsA[i].CRD_COLOR;
622 cb2 = pointsA[i1].CRD_COLOR;
623 cb3 = pointsB[ii].CRD_COLOR;
624 cb4 = pointsB[ii1].CRD_COLOR;
626 cb1 = z2cb(pointsA[i].z);
627 cb2 = z2cb(pointsA[i1].z);
628 cb3 = z2cb(pointsB[ii].z);
629 cb4 = z2cb(pointsB[ii1].z);
631 switch (pm3d.which_corner_color) {
632 case PM3D_WHICHCORNER_MEAN: avgC = (cb1 + cb2 + cb3 + cb4) * 0.25; break;
633 case PM3D_WHICHCORNER_GEOMEAN: avgC = geomean4(cb1, cb2, cb3, cb4); break;
634 case PM3D_WHICHCORNER_MEDIAN: avgC = median4(cb1, cb2, cb3, cb4); break;
635 case PM3D_WHICHCORNER_MIN: avgC = minimum4(cb1, cb2, cb3, cb4); break;
636 case PM3D_WHICHCORNER_MAX: avgC = maximum4(cb1, cb2, cb3, cb4); break;
637 case PM3D_WHICHCORNER_C1: avgC = cb1; break;
638 case PM3D_WHICHCORNER_C2: avgC = cb2; break;
639 case PM3D_WHICHCORNER_C3: avgC = cb3; break;
640 case PM3D_WHICHCORNER_C4: avgC = cb4; break;
641 default: int_error(NO_CARET, "cannot be here"); avgC = 0;
643 /* transform z value to gray, i.e. to interval [0,1] */
644 gray = cb2gray(avgC);
646 /* print the quadrangle with the given color */
647 printf("averageColor %g\tgray=%g\tM %g %g L %g %g L %g %g L %g %g\n",
650 pointsA[i].x, pointsA[i].y, pointsB[ii].x, pointsB[ii].y,
651 pointsB[ii1].x, pointsB[ii1].y, pointsA[i1].x, pointsA[i1].y);
654 if (pm3d.direction != PM3D_DEPTH)
656 #ifdef EXTENDED_COLOR_SPECS
661 corners[0].x = pointsA[i].x;
662 corners[0].y = pointsA[i].y;
663 corners[1].x = pointsB[ii].x;
664 corners[1].y = pointsB[ii].y;
665 corners[2].x = pointsB[ii1].x;
666 corners[2].y = pointsB[ii1].y;
667 corners[3].x = pointsA[i1].x;
668 corners[3].y = pointsA[i1].y;
670 if ( pm3d.interp_i > 1 || pm3d.interp_j > 1 || at_which_z == PM3D_AT_SURFACE) {
671 /* always supply the z value if
672 * EXTENDED_COLOR_SPECS is defined
674 corners[0].z = pointsA[i].z;
675 corners[1].z = pointsB[ii].z;
676 corners[2].z = pointsB[ii1].z;
677 corners[3].z = pointsA[i1].z;
678 if (color_from_column) {
679 corners[0].c = pointsA[i].CRD_COLOR;
680 corners[1].c = pointsB[ii].CRD_COLOR;
681 corners[2].c = pointsB[ii1].CRD_COLOR;
682 corners[3].c = pointsA[i1].CRD_COLOR;
685 #ifdef EXTENDED_COLOR_SPECS
686 if (supply_extended_color_specs) {
687 if (color_from_column) {
688 icorners[0].z = pointsA[i].CRD_COLOR;
689 icorners[1].z = pointsB[ii].CRD_COLOR;
690 icorners[2].z = pointsB[ii1].CRD_COLOR;
691 icorners[3].z = pointsA[i1].CRD_COLOR;
693 /* the target wants z and gray value */
694 icorners[0].z = pointsA[i].z;
695 icorners[1].z = pointsB[ii].z;
696 icorners[2].z = pointsB[ii1].z;
697 icorners[3].z = pointsA[i1].z;
699 for (i = 0; i < 4; i++) {
700 icorners[i].spec.gray =
701 cb2gray( color_from_column ? icorners[i].z : z2cb(icorners[i].z) );
704 if (pm3d.direction == PM3D_DEPTH) {
705 /* copy quadrangle */
706 quadrangle* qp = quadrangles + current_quadrangle;
707 memcpy(qp->corners, corners, 4 * sizeof (gpdPoint));
709 for (i = 0; i < 4; i++) {
710 qp->icorners[i].z = icorners[i].z;
711 qp->icorners[i].spec.gray = icorners[i].spec.gray;
713 current_quadrangle++;
716 filled_quadrangle(corners, icorners);
718 if (pm3d.interp_i > 1 || pm3d.interp_j > 1) {
719 /* Interpolation is enabled.
720 * interp_i is the # of points along scan lines
721 * interp_j is the # of points between scan lines
722 * Algorithm is to first sample i points along the scan lines
723 * defined by corners[3],corners[0] and corners[2],corners[1]. */
725 for (i1 = 0; i1 <= pm3d.interp_i; i1++) {
727 ((corners[3].x - corners[0].x) / pm3d.interp_i) * i1 + corners[0].x;
728 bl_point[i1][pm3d.interp_j].x =
729 ((corners[2].x - corners[1].x) / pm3d.interp_i) * i1 + corners[1].x;
731 ((corners[3].y - corners[0].y) / pm3d.interp_i) * i1 + corners[0].y;
732 bl_point[i1][pm3d.interp_j].y =
733 ((corners[2].y - corners[1].y) / pm3d.interp_i) * i1 + corners[1].y;
735 ((corners[3].z - corners[0].z) / pm3d.interp_i) * i1 + corners[0].z;
736 bl_point[i1][pm3d.interp_j].z =
737 ((corners[2].z - corners[1].z) / pm3d.interp_i) * i1 + corners[1].z;
738 if (color_from_column) {
740 ((corners[3].c - corners[0].c) / pm3d.interp_i) * i1 + corners[0].c;
741 bl_point[i1][pm3d.interp_j].c =
742 ((corners[2].c - corners[1].c) / pm3d.interp_i) * i1 + corners[1].c;
744 /* Next we sample j points between each of the new points
745 * created in the previous step (this samples between
746 * scan lines) in the same manner. */
747 for (j1 = 1; j1 < pm3d.interp_j; j1++) {
749 ((bl_point[i1][pm3d.interp_j].x - bl_point[i1][0].x) / pm3d.interp_j) *
750 j1 + bl_point[i1][0].x;
752 ((bl_point[i1][pm3d.interp_j].y - bl_point[i1][0].y) / pm3d.interp_j) *
753 j1 + bl_point[i1][0].y;
755 ((bl_point[i1][pm3d.interp_j].z - bl_point[i1][0].z) / pm3d.interp_j) *
756 j1 + bl_point[i1][0].z;
757 if (color_from_column)
759 ((bl_point[i1][pm3d.interp_j].c - bl_point[i1][0].c) / pm3d.interp_j) *
760 j1 + bl_point[i1][0].c;
763 /* Once all points are created, move them into an appropriate
764 * structure and call set_color on each to retrieve the
765 * correct color mapping for this new sub-sampled quadrangle. */
766 for (i1 = 0; i1 < pm3d.interp_i; i1++) {
767 for (j1 = 0; j1 < pm3d.interp_j; j1++) {
768 corners[0].x = bl_point[i1][j1].x;
769 corners[0].y = bl_point[i1][j1].y;
770 corners[0].z = bl_point[i1][j1].z;
771 corners[1].x = bl_point[i1+1][j1].x;
772 corners[1].y = bl_point[i1+1][j1].y;
773 corners[1].z = bl_point[i1+1][j1].z;
774 corners[2].x = bl_point[i1+1][j1+1].x;
775 corners[2].y = bl_point[i1+1][j1+1].y;
776 corners[2].z = bl_point[i1+1][j1+1].z;
777 corners[3].x = bl_point[i1][j1+1].x;
778 corners[3].y = bl_point[i1][j1+1].y;
779 corners[3].z = bl_point[i1][j1+1].z;
780 if (color_from_column) {
781 corners[0].c = bl_point[i1][j1].c;
782 corners[1].c = bl_point[i1+1][j1].c;
783 corners[2].c = bl_point[i1+1][j1+1].c;
784 corners[3].c = bl_point[i1][j1+1].c;
787 printf("(%g,%g),(%g,%g),(%g,%g),(%g,%g)\n",
788 corners[0].x, corners[0].y,
789 corners[1].x, corners[1].y,
790 corners[2].x, corners[2].y,
791 corners[3].x, corners[3].y);
793 /* If the colors are given separately, we already loaded them above */
794 if (!color_from_column) {
795 cb1 = z2cb(corners[0].z);
796 cb2 = z2cb(corners[1].z);
797 cb3 = z2cb(corners[2].z);
798 cb4 = z2cb(corners[3].z);
805 switch (pm3d.which_corner_color) {
806 case PM3D_WHICHCORNER_MEAN: avgC = (cb1 + cb2 + cb3 + cb4) * 0.25; break;
807 case PM3D_WHICHCORNER_GEOMEAN: avgC = geomean4(cb1, cb2, cb3, cb4); break;
808 case PM3D_WHICHCORNER_MEDIAN: avgC = median4(cb1, cb2, cb3, cb4); break;
809 case PM3D_WHICHCORNER_MIN: avgC = minimum4(cb1, cb2, cb3, cb4); break;
810 case PM3D_WHICHCORNER_MAX: avgC = maximum4(cb1, cb2, cb3, cb4); break;
811 case PM3D_WHICHCORNER_C1: avgC = cb1; break;
812 case PM3D_WHICHCORNER_C2: avgC = cb2; break;
813 case PM3D_WHICHCORNER_C3: avgC = cb3; break;
814 case PM3D_WHICHCORNER_C4: avgC = cb4; break;
815 default: int_error(NO_CARET, "cannot be here"); avgC = 0;
817 /* transform z value to gray, i.e. to interval [0,1] */
818 gray = cb2gray(avgC);
819 if (pm3d.direction != PM3D_DEPTH) {
821 filled_quadrangle(corners);
823 /* copy quadrangle */
824 quadrangle* qp = quadrangles + current_quadrangle;
825 memcpy(qp->corners, corners, 4 * sizeof (gpdPoint));
827 current_quadrangle++;
831 } else { /* thus (pm3d.interp_i == 1 && pm3d.interp_j == 1) */
832 if (pm3d.direction != PM3D_DEPTH) {
833 filled_quadrangle(corners);
835 /* copy quadrangle */
836 quadrangle* qp = quadrangles + current_quadrangle;
837 memcpy(qp->corners, corners, 4 * sizeof (gpdPoint));
839 current_quadrangle++;
841 } /* interpolate between points */
843 } /* loop quadrangles over points of two subsequent scans */
844 } /* loop over scans */
847 for (i1 = 0; i1 <= pm3d.interp_i; i1++)
851 /* free memory allocated by scan_array */
853 } /* end of pm3d splotting mode */
857 * Now the implementation of the filled color contour plot
860 filled_color_contour_plot(struct surface_points *this_plot, int contours_where)
863 struct gnuplot_contours *cntr;
865 /* just a shortcut */
866 TBOOLEAN color_from_column = this_plot->pm3d_color_from_column;
868 if (this_plot == NULL || this_plot->contours == NULL)
870 if (contours_where != CONTOUR_SRF && contours_where != CONTOUR_BASE)
873 /* return if the terminal does not support filled polygons */
874 if (!term->filled_polygon)
877 /* TODO: CHECK FOR NUMBER OF POINTS IN CONTOUR: IF TOO SMALL, THEN IGNORE! */
878 cntr = this_plot->contours;
880 printf("# Contour: points %i, z %g, label: %s\n", cntr->num_pts, cntr->coords[0].z, (cntr->label) ? cntr->label : "<no>");
881 if (cntr->isNewLevel) {
882 printf("\t...it isNewLevel\n");
883 /* contour split across chunks */
884 /* fprintf(gpoutfile, "\n# Contour %d, label: %s\n", number++, c->label); */
885 /* What is the color? */
886 /* get the z-coordinate */
887 /* transform contour z-coordinate value to gray, i.e. to interval [0,1] */
888 if (color_from_column)
889 gray = cb2gray(cntr->coords[0].CRD_COLOR);
891 gray = cb2gray( z2cb(cntr->coords[0].z) );
894 /* draw one countour */
895 if (contours_where == CONTOUR_SRF) /* at CONTOUR_SRF */
896 filled_polygon_3dcoords(cntr->num_pts, cntr->coords);
897 else /* at CONTOUR_BASE */
898 filled_polygon_3dcoords_zfixed(cntr->num_pts, cntr->coords, base_z);
902 } /* end of filled color contour plot splot mode */
906 * unset pm3d for the reset command
911 strcpy(pm3d.where, "s");
912 pm3d.flush = PM3D_FLUSH_BEGIN;
914 pm3d.direction = PM3D_SCANS_AUTOMATIC;
915 pm3d.clip = PM3D_CLIP_4IN;
916 pm3d.hidden3d_tag = 0;
919 #endif /* PM3D_HAVE_SOLID */
920 pm3d.implicit = PM3D_EXPLICIT;
921 pm3d.which_corner_color = PM3D_WHICHCORNER_MEAN;
928 * Draw (one) PM3D color surface.
931 pm3d_draw_one(struct surface_points *plot)
934 char *where = plot->pm3d_where[0] ? plot->pm3d_where : pm3d.where;
935 /* Draw either at 'where' option of the given surface or at pm3d.where
942 /* for pm3dCompress.awk */
943 if (gppsfile && (pm3d.direction != PM3D_DEPTH))
944 fputs("%pm3d_map_begin\n", gppsfile);
946 for (; where[i]; i++) {
947 pm3d_plot(plot, where[i]);
950 if (strchr(where, 'C') != NULL) {
951 /* !!!!! FILLED COLOR CONTOURS, *UNDOCUMENTED*
952 !!!!! LATER CHANGE TO STH LIKE
953 !!!!! (if_filled_contours_requested)
955 Currently filled color contours do not work because gnuplot generates
956 open contour lines, i.e. not closed on the graph boundary.
958 if (draw_contour & CONTOUR_SRF)
959 filled_color_contour_plot(plot, CONTOUR_SRF);
960 if (draw_contour & CONTOUR_BASE)
961 filled_color_contour_plot(plot, CONTOUR_BASE);
964 /* for pm3dCompress.awk */
965 if (gppsfile && (pm3d.direction != PM3D_DEPTH))
966 fputs("%pm3d_map_end\n", gppsfile);
970 /* Display an error message for the routine get_pm3d_at_option() below.
974 pm3d_option_at_error()
976 int_error(c_token, "\
977 parameter to `pm3d at` requires combination of up to 6 characters b,s,t\n\
978 \t(drawing at bottom, surface, top)");
982 /* Read the option for 'pm3d at' command.
983 * Used by 'set pm3d at ...' or by 'splot ... with pm3d at ...'.
984 * If no option given, then returns empty string, otherwise copied there.
985 * The string is unchanged on error, and 1 is returned.
986 * On success, 0 is returned.
989 get_pm3d_at_option(char *pm3d_where)
993 if (END_OF_COMMAND || token[c_token].length >= sizeof(pm3d.where)) {
994 pm3d_option_at_error();
997 memcpy(pm3d_where, gp_input_line + token[c_token].start_index,
998 token[c_token].length);
999 pm3d_where[ token[c_token].length ] = 0;
1000 /* verify the parameter */
1001 for (c = pm3d_where; *c; c++) {
1002 if (*c != 'C') /* !!!!! CONTOURS, UNDOCUMENTED, THIS LINE IS TEMPORARILY HERE !!!!! */
1003 if (*c != PM3D_AT_BASE && *c != PM3D_AT_TOP && *c != PM3D_AT_SURFACE) {
1004 pm3d_option_at_error();
1012 /* Set flag plot_has_palette to TRUE if there is any element on the graph
1013 * which requires palette of continuous colors.
1016 set_plot_with_palette(int plot_num, int plot_mode)
1018 struct surface_points *this_3dplot = first_3dplot;
1019 struct curve_points *this_2dplot = first_plot;
1021 struct text_label *this_label = first_label;
1023 plot_has_palette = TRUE;
1024 /* Is pm3d switched on globally? */
1025 if (pm3d.implicit == PM3D_IMPLICIT)
1028 /* Check 2D plots */
1029 if (plot_mode == MODE_PLOT) {
1030 while (this_2dplot) {
1032 if (this_2dplot->plot_style == IMAGE)
1035 if (this_2dplot->lp_properties.use_palette
1036 && this_2dplot->lp_properties.pm3d_color.type > TC_RGB)
1038 if (this_2dplot->labels
1039 && this_2dplot->labels->textcolor.type >= TC_CB)
1041 this_2dplot = this_2dplot->next;
1045 /* Check 3D plots */
1046 if (plot_mode == MODE_SPLOT) {
1047 /* Any surface 'with pm3d', 'with image' or 'with line|dot palette'? */
1048 while (surface++ < plot_num) {
1049 if (this_3dplot->plot_style == PM3DSURFACE)
1052 if (this_3dplot->plot_style == IMAGE)
1055 if (this_3dplot->lp_properties.use_palette) {
1056 int type = this_3dplot->lp_properties.pm3d_color.type;
1057 if (type == TC_LT || type == TC_LINESTYLE || type == TC_RGB)
1058 /* don't return yet */
1061 /* TC_DEFAULT: splot x with line|lp|dot palette */
1064 #ifdef EAM_DATASTRINGS
1065 if (this_3dplot->labels &&
1066 this_3dplot->labels->textcolor.type >= TC_CB)
1069 this_3dplot = this_3dplot->next_sp;
1073 /* Any label with 'textcolor palette'? */
1074 #define TC_USES_PALETTE(tctype) (tctype==TC_Z) || (tctype==TC_CB) || (tctype==TC_FRAC)
1075 for (; this_label != NULL; this_label = this_label->next) {
1076 if (TC_USES_PALETTE(this_label->textcolor.type))
1079 /* Any of title, xlabel, ylabel, zlabel, ... with 'textcolor palette'? */
1080 if (TC_USES_PALETTE(title.textcolor.type)) return;
1081 if (TC_USES_PALETTE(axis_array[FIRST_X_AXIS].label.textcolor.type)) return;
1082 if (TC_USES_PALETTE(axis_array[FIRST_Y_AXIS].label.textcolor.type)) return;
1083 if (TC_USES_PALETTE(axis_array[SECOND_X_AXIS].label.textcolor.type)) return;
1084 if (TC_USES_PALETTE(axis_array[SECOND_Y_AXIS].label.textcolor.type)) return;
1085 if (plot_mode == MODE_SPLOT)
1086 if (TC_USES_PALETTE(axis_array[FIRST_Z_AXIS].label.textcolor.type)) return;
1087 if (TC_USES_PALETTE(axis_array[COLOR_AXIS].label.textcolor.type)) return;
1088 #undef TC_USES_PALETTE
1090 /* Palette with continuous colors is not used. */
1091 plot_has_palette = FALSE; /* otherwise it stays TRUE */
1095 is_plot_with_palette()
1097 return plot_has_palette;
1101 is_plot_with_colorbox()
1103 return plot_has_palette && (color_box.where != SMCOLOR_BOX_NO);