1 /* GNUPLOT - pm3d.c */
2
3 /*[
4 *
5 * Petr Mikulik, since December 1998
6 * Copyright: open source as much as possible
7 *
8 * What is here: global variables and routines for the pm3d splotting mode.
9 *
10 ]*/
11
12 #ifdef HAVE_CONFIG_H
13 # include "config.h"
14 #endif
15 #include "pm3d.h"
16 #include "alloc.h"
17 #include "getcolor.h" /* for rgb1_from_gray */
18 #include "graphics.h"
19 #include "hidden3d.h" /* p_vertex & map3d_xyz() */
20 #include "plot2d.h"
21 #include "plot3d.h"
22 #include "setshow.h" /* for surface_rot_z */
23 #include "term_api.h" /* for lp_use_properties() */
24 #include "command.h" /* for c_token */
25 #include "voxelgrid.h" /* for isosurface_options */
26
27 #include <stdlib.h> /* qsort() */
28
29 /* Used to initialize `set pm3d border` */
30 struct lp_style_type default_pm3d_border = DEFAULT_LP_STYLE_TYPE;
31
32 /* Set by plot styles that use pm3d quadrangles even in non-pm3d mode */
33 TBOOLEAN track_pm3d_quadrangles;
34
35 /*
36 Global options for pm3d algorithm (to be accessed by set / show).
37 */
38 pm3d_struct pm3d; /* Initialized via init_session->reset_command->reset_pm3d */
39
40 lighting_model pm3d_shade;
41
42 typedef struct {
43 double gray;
44 double z; /* z value after rotation to graph coordinates for depth sorting */
45 union {
46 gpdPoint corners[4]; /* The normal case. vertices stored right here */
47 int array_index; /* Stored elsewhere if there are more than 4 */
48 } vertex;
49 union { /* Only used by depthorder processing */
50 t_colorspec *colorspec;
51 unsigned int rgb_color;
52 } qcolor;
53 short fillstyle; /* from plot->fill_properties */
54 short type; /* QUAD_TYPE_NORMAL or QUAD_TYPE_4SIDES etc */
55 } quadrangle;
56
57 #define QUAD_TYPE_NORMAL 0
58 #define QUAD_TYPE_TRIANGLE 3
59 #define QUAD_TYPE_4SIDES 4
60 #define QUAD_TYPE_LARGEPOLYGON 5
61
62 #define PM3D_USE_COLORSPEC_INSTEAD_OF_GRAY -12345
63 #define PM3D_USE_RGB_COLOR_INSTEAD_OF_GRAY -12346
64 #define PM3D_USE_BACKGROUND_INSTEAD_OF_GRAY -12347
65
66 static int allocated_quadrangles = 0;
67 static int current_quadrangle = 0;
68 static quadrangle* quadrangles = (quadrangle*)0;
69
70 static gpdPoint *polygonlist = NULL; /* holds polygons with >4 vertices */
71 static int next_polygon = 0; /* index of next slot in the list */
72 static int current_polygon = 0; /* index of the current polygon */
73 static int polygonlistsize = 0;
74
75 static int pm3d_plot_at = 0; /* flag so that top/base polygons are not clipped against z */
76
77 /* Internal prototypes for this module */
78 static TBOOLEAN plot_has_palette;
79 static double geomean4(double, double, double, double);
80 static double harmean4(double, double, double, double);
81 static double median4(double, double, double, double);
82 static double rms4(double, double, double, double);
83 static void pm3d_plot(struct surface_points *, int);
84 static void pm3d_option_at_error(void);
85 static void pm3d_rearrange_part(struct iso_curve *, const int, struct iso_curve ***, int *);
86 static int apply_lighting_model( struct coordinate *, struct coordinate *, struct coordinate *, struct coordinate *, double gray );
87 static void illuminate_one_quadrangle( quadrangle *q );
88
89 static void filled_polygon(gpdPoint *corners, int fillstyle, int nv);
90 static int clip_filled_polygon( gpdPoint *inpts, gpdPoint *outpts, int nv );
91
92 static TBOOLEAN color_from_rgbvar = FALSE;
93 static double light[3];
94
95 static gpdPoint *get_polygon(int size);
96 static void free_polygonlist(void);
97
98 /*
99 * Utility routines.
100 */
101
102 /* Geometrical mean = pow( prod(x_i > 0) x_i, 1/N )
103 * In order to facilitate plotting surfaces with all color coords negative,
104 * All 4 corners positive - return positive geometric mean
105 * All 4 corners negative - return negative geometric mean
106 * otherwise return 0
107 * EAM Oct 2012: This is a CHANGE
108 */
109 static double
geomean4(double x1,double x2,double x3,double x4)110 geomean4 (double x1, double x2, double x3, double x4)
111 {
112 int neg = (x1 < 0) + (x2 < 0) + (x3 < 0) + (x4 < 0);
113 double product = x1 * x2 * x3 * x4;
114
115 if (product == 0) return 0;
116 if (neg == 1 || neg == 2 || neg == 3) return 0;
117
118 product = sqrt(sqrt(fabs(product)));
119 return (neg == 0) ? product : -product;
120 }
121
122 static double
harmean4(double x1,double x2,double x3,double x4)123 harmean4 (double x1, double x2, double x3, double x4)
124 {
125 if (x1 <= 0 || x2 <= 0 || x3 <= 0 || x4 <= 0)
126 return not_a_number();
127 else
128 return 4 / ((1/x1) + (1/x2) + (1/x3) + (1/x4));
129 }
130
131
132 /* Median: sort values, and then: for N odd, it is the middle value; for N even,
133 * it is mean of the two middle values.
134 */
135 static double
median4(double x1,double x2,double x3,double x4)136 median4 (double x1, double x2, double x3, double x4)
137 {
138 double tmp;
139 /* sort them: x1 < x2 and x3 < x4 */
140 if (x1 > x2) { tmp = x2; x2 = x1; x1 = tmp; }
141 if (x3 > x4) { tmp = x3; x3 = x4; x4 = tmp; }
142 /* sum middle numbers */
143 tmp = (x1 < x3) ? x3 : x1;
144 tmp += (x2 < x4) ? x2 : x4;
145 return tmp * 0.5;
146 }
147
148
149 /* Minimum of 4 numbers.
150 */
151 static double
minimum4(double x1,double x2,double x3,double x4)152 minimum4 (double x1, double x2, double x3, double x4)
153 {
154 x1 = GPMIN(x1, x2);
155 x3 = GPMIN(x3, x4);
156 return GPMIN(x1, x3);
157 }
158
159
160 /* Maximum of 4 numbers.
161 */
162 static double
maximum4(double x1,double x2,double x3,double x4)163 maximum4 (double x1, double x2, double x3, double x4)
164 {
165 x1 = GPMAX(x1, x2);
166 x3 = GPMAX(x3, x4);
167 return GPMAX(x1, x3);
168 }
169
170 /* The root mean square of the 4 numbers */
171 static double
rms4(double x1,double x2,double x3,double x4)172 rms4 (double x1, double x2, double x3, double x4)
173 {
174 return 0.5*sqrt(x1*x1 + x2*x2 + x3*x3 + x4*x4);
175 }
176
177 /*
178 * Now the routines which are really just those for pm3d.c
179 */
180
181 /*
182 * Rescale cb (color) value into the interval of grays [0,1], taking care
183 * of palette being positive or negative.
184 * Note that it is OK for logarithmic cb-axis too.
185 */
186 double
cb2gray(double cb)187 cb2gray(double cb)
188 {
189 AXIS *cbaxis = &CB_AXIS;
190
191 if (cb <= cbaxis->min)
192 return (sm_palette.positive == SMPAL_POSITIVE) ? 0 : 1;
193 if (cb >= cbaxis->max)
194 return (sm_palette.positive == SMPAL_POSITIVE) ? 1 : 0;
195
196 if (nonlinear(cbaxis)) {
197 cbaxis = cbaxis->linked_to_primary;
198 cb = eval_link_function(cbaxis, cb);
199 }
200
201 cb = (cb - cbaxis->min) / (cbaxis->max - cbaxis->min);
202 return (sm_palette.positive == SMPAL_POSITIVE) ? cb : 1-cb;
203 }
204
205
206 /*
207 * Rearrange...
208 */
209 static void
pm3d_rearrange_part(struct iso_curve * src,const int len,struct iso_curve *** dest,int * invert)210 pm3d_rearrange_part(
211 struct iso_curve *src,
212 const int len,
213 struct iso_curve ***dest,
214 int *invert)
215 {
216 struct iso_curve *scanA;
217 struct iso_curve *scanB;
218 struct iso_curve **scan_array;
219 int i, scan;
220 int invert_order = 0;
221
222 /* loop over scans in one surface
223 Scans are linked from this_plot->iso_crvs in the opposite order than
224 they are in the datafile.
225 Therefore it is necessary to make vector scan_array of iso_curves.
226 Scans are sorted in scan_array according to pm3d.direction (this can
227 be PM3D_SCANS_FORWARD or PM3D_SCANS_BACKWARD).
228 */
229 scan_array = *dest = gp_alloc(len * sizeof(scanA), "pm3d scan array");
230
231 if (pm3d.direction == PM3D_SCANS_AUTOMATIC) {
232 int cnt;
233 int len2 = len;
234 TBOOLEAN exit_outer_loop = 0;
235
236 for (scanA = src; scanA && 0 == exit_outer_loop; scanA = scanA->next, len2--) {
237
238 int from, i;
239 vertex vA, vA2;
240
241 if ((cnt = scanA->p_count - 1) <= 0)
242 continue;
243
244 /* ordering within one scan */
245 for (from=0; from<=cnt; from++) /* find 1st non-undefined point */
246 if (scanA->points[from].type != UNDEFINED) {
247 map3d_xyz(scanA->points[from].x, scanA->points[from].y, 0, &vA);
248 break;
249 }
250 for (i=cnt; i>from; i--) /* find the last non-undefined point */
251 if (scanA->points[i].type != UNDEFINED) {
252 map3d_xyz(scanA->points[i].x, scanA->points[i].y, 0, &vA2);
253 break;
254 }
255
256 if (i - from > cnt * 0.1)
257 /* it is completely arbitrary to request at least
258 * 10% valid samples in this scan. (joze Jun-05-2002) */
259 *invert = (vA2.z > vA.z) ? 0 : 1;
260 else
261 continue; /* all points were undefined, so check next scan */
262
263
264 /* check the z ordering between scans
265 * Find last scan. If this scan has all points undefined,
266 * find last but one scan, an so on. */
267
268 for (; len2 >= 3 && !exit_outer_loop; len2--) {
269 for (scanB = scanA-> next, i = len2 - 2; i && scanB; i--)
270 scanB = scanB->next; /* skip over to last scan */
271 if (scanB && scanB->p_count) {
272 vertex vB;
273 for (i = from /* we compare vA.z with vB.z */; i<scanB->p_count; i++) {
274 /* find 1st non-undefined point */
275 if (scanB->points[i].type != UNDEFINED) {
276 map3d_xyz(scanB->points[i].x, scanB->points[i].y, 0, &vB);
277 invert_order = (vB.z > vA.z) ? 0 : 1;
278 exit_outer_loop = 1;
279 break;
280 }
281 }
282 }
283 }
284 }
285 }
286
287 FPRINTF((stderr, "(pm3d_rearrange_part) invert = %d\n", *invert));
288 FPRINTF((stderr, "(pm3d_rearrange_part) invert_order = %d\n", invert_order));
289
290 for (scanA = src, scan = len - 1, i = 0; scan >= 0; --scan, i++) {
291 if (pm3d.direction == PM3D_SCANS_AUTOMATIC) {
292 switch (invert_order) {
293 case 1:
294 scan_array[scan] = scanA;
295 break;
296 case 0:
297 default:
298 scan_array[i] = scanA;
299 break;
300 }
301 } else if (pm3d.direction == PM3D_SCANS_FORWARD)
302 scan_array[scan] = scanA;
303 else /* PM3D_SCANS_BACKWARD: i counts scans */
304 scan_array[i] = scanA;
305 scanA = scanA->next;
306 }
307 }
308
309
310 /*
311 * Rearrange scan array
312 *
313 * Allocates *first_ptr (and eventually *second_ptr)
314 * which must be freed by the caller
315 */
316 void
pm3d_rearrange_scan_array(struct surface_points * this_plot,struct iso_curve *** first_ptr,int * first_n,int * first_invert,struct iso_curve *** second_ptr,int * second_n,int * second_invert)317 pm3d_rearrange_scan_array(
318 struct surface_points *this_plot,
319 struct iso_curve ***first_ptr,
320 int *first_n, int *first_invert,
321 struct iso_curve ***second_ptr,
322 int *second_n, int *second_invert)
323 {
324 if (first_ptr) {
325 pm3d_rearrange_part(this_plot->iso_crvs, this_plot->num_iso_read, first_ptr, first_invert);
326 *first_n = this_plot->num_iso_read;
327 }
328 if (second_ptr) {
329 struct iso_curve *icrvs = this_plot->iso_crvs;
330 struct iso_curve *icrvs2;
331 int i;
332 /* advance until second part */
333 for (i = 0; i < this_plot->num_iso_read; i++)
334 icrvs = icrvs->next;
335 /* count the number of scans of second part */
336 for (i = 0, icrvs2 = icrvs; icrvs2; icrvs2 = icrvs2->next)
337 i++;
338 if (i > 0) {
339 *second_n = i;
340 pm3d_rearrange_part(icrvs, i, second_ptr, second_invert);
341 } else {
342 *second_ptr = (struct iso_curve **) 0;
343 }
344 }
345 }
346
compare_quadrangles(const void * v1,const void * v2)347 static int compare_quadrangles(const void* v1, const void* v2)
348 {
349 const quadrangle* q1 = (const quadrangle*)v1;
350 const quadrangle* q2 = (const quadrangle*)v2;
351
352 if (q1->z > q2->z)
353 return 1;
354 else if (q1->z < q2->z)
355 return -1;
356 else
357 return 0;
358 }
359
pm3d_depth_queue_clear(void)360 void pm3d_depth_queue_clear(void)
361 {
362 free(quadrangles);
363 quadrangles = NULL;
364 allocated_quadrangles = 0;
365 current_quadrangle = 0;
366 }
367
pm3d_depth_queue_flush(void)368 void pm3d_depth_queue_flush(void)
369 {
370 if (pm3d.direction != PM3D_DEPTH && !track_pm3d_quadrangles)
371 return;
372
373 term->layer(TERM_LAYER_BEGIN_PM3D_FLUSH);
374
375 if (current_quadrangle > 0 && quadrangles) {
376
377 quadrangle* qp;
378 quadrangle* qe;
379 gpdPoint* gpdPtr;
380 vertex out;
381 double zbase = 0;
382 int nv, i;
383
384 if (pm3d.base_sort)
385 cliptorange(zbase, Z_AXIS.min, Z_AXIS.max);
386
387 for (qp = quadrangles, qe = quadrangles + current_quadrangle; qp != qe; qp++) {
388 double z = 0;
389 double zmean = 0;
390
391 if (qp->type == QUAD_TYPE_LARGEPOLYGON) {
392 gpdPtr = &polygonlist[qp->vertex.array_index];
393 nv = gpdPtr[2].c;
394 } else {
395 gpdPtr = qp->vertex.corners;
396 nv = 4;
397 }
398
399
400 for (i = 0; i < nv; i++, gpdPtr++) {
401 /* 3D boxes want to be sorted on z of the base, not the mean z */
402 if (pm3d.base_sort)
403 map3d_xyz(gpdPtr->x, gpdPtr->y, zbase, &out);
404 else
405 map3d_xyz(gpdPtr->x, gpdPtr->y, gpdPtr->z, &out);
406 zmean += out.z;
407 if (i == 0 || out.z > z)
408 z = out.z;
409 }
410
411 if (pm3d.zmean_sort)
412 qp->z = zmean / nv;
413 else
414 qp->z = z;
415 }
416
417 qsort(quadrangles, current_quadrangle, sizeof (quadrangle), compare_quadrangles);
418
419 for (qp = quadrangles, qe = quadrangles + current_quadrangle; qp != qe; qp++) {
420
421 /* set the color */
422 if (qp->gray == PM3D_USE_COLORSPEC_INSTEAD_OF_GRAY)
423 apply_pm3dcolor(qp->qcolor.colorspec);
424 else if (qp->gray == PM3D_USE_BACKGROUND_INSTEAD_OF_GRAY)
425 term->linetype(LT_BACKGROUND);
426 else if (qp->gray == PM3D_USE_RGB_COLOR_INSTEAD_OF_GRAY)
427 set_rgbcolor_var(qp->qcolor.rgb_color);
428 else if (pm3d_shade.strength > 0)
429 set_rgbcolor_var((unsigned int)qp->gray);
430 else
431 set_color(qp->gray);
432
433 if (qp->type == QUAD_TYPE_LARGEPOLYGON) {
434 gpdPoint *vertices = &polygonlist[qp->vertex.array_index];
435 int nv = vertices[2].c;
436 filled_polygon(vertices, qp->fillstyle, nv);
437 } else {
438 filled_polygon(qp->vertex.corners, qp->fillstyle, 4);
439 }
440 }
441 }
442
443 pm3d_depth_queue_clear();
444 free_polygonlist();
445
446 term->layer(TERM_LAYER_END_PM3D_FLUSH);
447 }
448
449 /*
450 * Now the implementation of the pm3d (s)plotting mode
451 */
452 static void
pm3d_plot(struct surface_points * this_plot,int at_which_z)453 pm3d_plot(struct surface_points *this_plot, int at_which_z)
454 {
455 int j, i, i1, ii, ii1, from, scan, up_to, up_to_minus, invert = 0;
456 int go_over_pts, max_pts;
457 int are_ftriangles, ftriangles_low_pt = -999, ftriangles_high_pt = -999;
458 struct iso_curve *scanA, *scanB;
459 struct coordinate *pointsA, *pointsB;
460 struct iso_curve **scan_array;
461 int scan_array_n;
462 double avgC, gray = 0;
463 double cb1, cb2, cb3, cb4;
464 gpdPoint corners[4];
465 int interp_i, interp_j;
466 gpdPoint **bl_point = NULL; /* used for bilinear interpolation */
467 TBOOLEAN color_from_column = FALSE;
468 TBOOLEAN color_from_fillcolor = FALSE;
469 int plot_fillstyle;
470
471 /* should never happen */
472 if (this_plot == NULL)
473 return;
474
475 /* just a shortcut */
476 plot_fillstyle = style_from_fill(&this_plot->fill_properties);
477 color_from_column = this_plot->pm3d_color_from_column;
478 color_from_rgbvar = FALSE;
479
480 if (this_plot->lp_properties.pm3d_color.type == TC_RGB
481 && this_plot->lp_properties.pm3d_color.value == -1)
482 color_from_rgbvar = TRUE;
483
484 if (this_plot->fill_properties.border_color.type == TC_RGB
485 || this_plot->fill_properties.border_color.type == TC_LINESTYLE) {
486 color_from_rgbvar = TRUE;
487 color_from_fillcolor = TRUE;
488 }
489
490 /* Apply and save the user-requested line properties */
491 term_apply_lp_properties(&this_plot->lp_properties);
492
493 if (at_which_z != PM3D_AT_BASE && at_which_z != PM3D_AT_TOP && at_which_z != PM3D_AT_SURFACE)
494 return;
495 else
496 pm3d_plot_at = at_which_z; /* clip_filled_polygon() will check this */
497
498 /* return if the terminal does not support filled polygons */
499 if (!term->filled_polygon)
500 return;
501
502 /* for pm3dCompress.awk and pm3dConvertToImage.awk */
503 if (pm3d.direction != PM3D_DEPTH)
504 term->layer(TERM_LAYER_BEGIN_PM3D_MAP);
505
506 switch (at_which_z) {
507 case PM3D_AT_BASE:
508 corners[0].z = corners[1].z = corners[2].z = corners[3].z = base_z;
509 break;
510 case PM3D_AT_TOP:
511 corners[0].z = corners[1].z = corners[2].z = corners[3].z = ceiling_z;
512 break;
513 /* the 3rd possibility is surface, PM3D_AT_SURFACE, coded below */
514 }
515
516 scanA = this_plot->iso_crvs;
517
518 pm3d_rearrange_scan_array(this_plot, &scan_array, &scan_array_n, &invert, (struct iso_curve ***) 0, (int *) 0, (int *) 0);
519
520 interp_i = pm3d.interp_i;
521 interp_j = pm3d.interp_j;
522
523 if (interp_i <= 0 || interp_j <= 0) {
524 /* Number of interpolations will be determined from desired number of points.
525 Search for number of scans and maximal number of points in a scan for points
526 which will be plotted (INRANGE). Then set interp_i,j so that number of points
527 will be a bit larger than |interp_i,j|.
528 If (interp_i,j==0) => set this number of points according to DEFAULT_OPTIMAL_NB_POINTS.
529 Ideally this should be comparable to the resolution of the output device, which
530 can hardly by done at this high level instead of the driver level.
531 */
532 #define DEFAULT_OPTIMAL_NB_POINTS 200
533 int max_scan_pts = 0;
534 int max_scans = 0;
535 int pts;
536 for (scan = 0; scan < this_plot->num_iso_read - 1; scan++) {
537 scanA = scan_array[scan];
538 pointsA = scanA->points;
539 pts = 0;
540 for (j=0; j<scanA->p_count; j++)
541 if (pointsA[j].type == INRANGE) pts++;
542 if (pts > 0) {
543 max_scan_pts = GPMAX(max_scan_pts, pts);
544 max_scans++;
545 }
546 }
547
548 if (max_scan_pts == 0 || max_scans == 0)
549 int_error(NO_CARET, "all scans empty");
550
551 if (interp_i <= 0) {
552 ii = (interp_i == 0) ? DEFAULT_OPTIMAL_NB_POINTS : -interp_i;
553 interp_i = floor(ii / max_scan_pts) + 1;
554 }
555 if (interp_j <= 0) {
556 ii = (interp_j == 0) ? DEFAULT_OPTIMAL_NB_POINTS : -interp_j;
557 interp_j = floor(ii / max_scans) + 1;
558 }
559 #if 0
560 fprintf(stderr, "pm3d.interp_i=%i\t pm3d.interp_j=%i\n", pm3d.interp_i, pm3d.interp_j);
561 fprintf(stderr, "INRANGE: max_scans=%i max_scan_pts=%i\n", max_scans, max_scan_pts);
562 fprintf(stderr, "setting interp_i=%i\t interp_j=%i => there will be %i and %i points\n",
563 interp_i, interp_j, interp_i*max_scan_pts, interp_j*max_scans);
564 #endif
565 }
566
567 if (pm3d.direction == PM3D_DEPTH) {
568 int needed_quadrangles = 0;
569
570 for (scan = 0; scan < this_plot->num_iso_read - 1; scan++) {
571
572 scanA = scan_array[scan];
573 scanB = scan_array[scan + 1];
574
575 are_ftriangles = pm3d.ftriangles && (scanA->p_count != scanB->p_count);
576 if (!are_ftriangles)
577 needed_quadrangles += GPMIN(scanA->p_count, scanB->p_count) - 1;
578 else {
579 needed_quadrangles += GPMAX(scanA->p_count, scanB->p_count) - 1;
580 }
581 }
582 needed_quadrangles *= (interp_i > 1) ? interp_i : 1;
583 needed_quadrangles *= (interp_j > 1) ? interp_j : 1;
584
585 if (needed_quadrangles > 0) {
586 while (current_quadrangle + needed_quadrangles >= allocated_quadrangles) {
587 FPRINTF((stderr, "allocated_quadrangles = %d current = %d needed = %d\n",
588 allocated_quadrangles, current_quadrangle, needed_quadrangles));
589 allocated_quadrangles = needed_quadrangles + 2*allocated_quadrangles;
590 quadrangles = (quadrangle*)gp_realloc(quadrangles,
591 allocated_quadrangles * sizeof (quadrangle),
592 "pm3d_plot->quadrangles");
593 }
594 }
595 }
596 /* pm3d_rearrange_scan_array(this_plot, (struct iso_curve***)0, (int*)0, &scan_array, &invert); */
597
598 #if 0
599 /* debugging: print scan_array */
600 for (scan = 0; scan < this_plot->num_iso_read; scan++) {
601 printf("**** SCAN=%d points=%d\n", scan, scan_array[scan]->p_count);
602 }
603 #endif
604
605 #if 0
606 /* debugging: this loop prints properties of all scans */
607 for (scan = 0; scan < this_plot->num_iso_read; scan++) {
608 struct coordinate *points;
609 scanA = scan_array[scan];
610 printf("\n#IsoCurve = scan nb %d, %d points\n#x y z type(in,out,undef)\n", scan, scanA->p_count);
611 for (i = 0, points = scanA->points; i < scanA->p_count; i++) {
612 printf("%g %g %g %c\n",
613 points[i].x, points[i].y, points[i].z, points[i].type == INRANGE ? 'i' : points[i].type == OUTRANGE ? 'o' : 'u');
614 /* Note: INRANGE, OUTRANGE, UNDEFINED */
615 }
616 }
617 printf("\n");
618 #endif
619
620 /*
621 * if bilinear interpolation is enabled, allocate memory for the
622 * interpolated points here
623 */
624 if (interp_i > 1 || interp_j > 1) {
625 bl_point = (gpdPoint **)gp_alloc(sizeof(gpdPoint*) * (interp_i+1), "bl-interp along scan");
626 for (i1 = 0; i1 <= interp_i; i1++)
627 bl_point[i1] = (gpdPoint *)gp_alloc(sizeof(gpdPoint) * (interp_j+1), "bl-interp between scan");
628 }
629
630 /*
631 * this loop does the pm3d draw of joining two curves
632 *
633 * How the loop below works:
634 * - scanB = scan last read; scanA = the previous one
635 * - link the scan from A to B, then move B to A, then read B, then draw
636 */
637 for (scan = 0; scan < this_plot->num_iso_read - 1; scan++) {
638 scanA = scan_array[scan];
639 scanB = scan_array[scan + 1];
640
641 FPRINTF((stderr,"\n#IsoCurveA = scan nb %d has %d points ScanB has %d points\n",
642 scan, scanA->p_count, scanB->p_count));
643
644 pointsA = scanA->points;
645 pointsB = scanB->points;
646 /* if the number of points in both scans is not the same, then the
647 * starting index (offset) of scan B according to the flushing setting
648 * has to be determined
649 */
650 from = 0; /* default is pm3d.flush==PM3D_FLUSH_BEGIN */
651 if (pm3d.flush == PM3D_FLUSH_END)
652 from = abs(scanA->p_count - scanB->p_count);
653 else if (pm3d.flush == PM3D_FLUSH_CENTER)
654 from = abs(scanA->p_count - scanB->p_count) / 2;
655 /* find the minimal number of points in both scans */
656 up_to = GPMIN(scanA->p_count, scanB->p_count) - 1;
657 up_to_minus = up_to - 1; /* calculate only once */
658 are_ftriangles = pm3d.ftriangles && (scanA->p_count != scanB->p_count);
659 if (!are_ftriangles)
660 go_over_pts = up_to;
661 else {
662 max_pts = GPMAX(scanA->p_count, scanB->p_count);
663 go_over_pts = max_pts - 1;
664 /* the j-subrange of quadrangles; in the remainder of the interval
665 * [0..up_to] the flushing triangles are to be drawn */
666 ftriangles_low_pt = from;
667 ftriangles_high_pt = from + up_to_minus;
668 }
669 /* Go over
670 * - the minimal number of points from both scans, if only quadrangles.
671 * - the maximal number of points from both scans if flush triangles
672 * (the missing points in the scan of lower nb of points will be
673 * duplicated from the begin/end points).
674 *
675 * Notice: if it would be once necessary to go over points in `backward'
676 * direction, then the loop body below would require to replace the data
677 * point indices `i' by `up_to-i' and `i+1' by `up_to-i-1'.
678 */
679 for (j = 0; j < go_over_pts; j++) {
680 /* Now i be the index of the scan with smaller number of points,
681 * ii of the scan with larger number of points. */
682 if (are_ftriangles && (j < ftriangles_low_pt || j > ftriangles_high_pt)) {
683 i = (j <= ftriangles_low_pt) ? 0 : ftriangles_high_pt-from+1;
684 ii = j;
685 i1 = i;
686 ii1 = ii + 1;
687 } else {
688 int jj = are_ftriangles ? j - from : j;
689 i = jj;
690 if (PM3D_SCANS_AUTOMATIC == pm3d.direction && invert)
691 i = up_to_minus - jj;
692 ii = i + from;
693 i1 = i + 1;
694 ii1 = ii + 1;
695 }
696 /* From here, i is index to scan A, ii to scan B */
697 if (scanA->p_count > scanB->p_count) {
698 int itmp = i;
699 i = ii;
700 ii = itmp;
701 itmp = i1;
702 i1 = ii1;
703 ii1 = itmp;
704 }
705
706 FPRINTF((stderr,"j=%i: i=%i i1=%i [%i] ii=%i ii1=%i [%i]\n",
707 j,i,i1,scanA->p_count,ii,ii1,scanB->p_count));
708
709 /* choose the clipping method */
710 if (pm3d.clip == PM3D_CLIP_4IN) {
711 /* (1) all 4 points of the quadrangle must be in x and y range */
712 if (!(pointsA[i].type == INRANGE && pointsA[i1].type == INRANGE &&
713 pointsB[ii].type == INRANGE && pointsB[ii1].type == INRANGE))
714 continue;
715 } else {
716 /* (2) all 4 points of the quadrangle must be defined */
717 if (pointsA[i].type == UNDEFINED || pointsA[i1].type == UNDEFINED ||
718 pointsB[ii].type == UNDEFINED || pointsB[ii1].type == UNDEFINED)
719 continue;
720 /* and at least 1 point of the quadrangle must be in x and y range */
721 if (pointsA[i].type == OUTRANGE && pointsA[i1].type == OUTRANGE &&
722 pointsB[ii].type == OUTRANGE && pointsB[ii1].type == OUTRANGE)
723 continue;
724 }
725
726 if ((interp_i <= 1 && interp_j <= 1) || pm3d.direction == PM3D_DEPTH) {
727 {
728 /* Get the gray as the average of the corner z- or gray-positions
729 (note: log scale is already included). The average is calculated here
730 if there is no interpolation (including the "pm3d depthorder" option),
731 otherwise it is done for each interpolated quadrangle later.
732 */
733 if (color_from_column) {
734 /* color is set in plot3d.c:get_3ddata() */
735 cb1 = pointsA[i].CRD_COLOR;
736 cb2 = pointsA[i1].CRD_COLOR;
737 cb3 = pointsB[ii].CRD_COLOR;
738 cb4 = pointsB[ii1].CRD_COLOR;
739 } else if (color_from_fillcolor) {
740 /* color is set by "fc <rgbvalue>" */
741 cb1 = cb2 = cb3 = cb4 = this_plot->fill_properties.border_color.lt;
742 /* EXPERIMENTAL
743 * pm3d fc linestyle N generates
744 * top/bottom color difference as with hidden3d
745 */
746 if (this_plot->fill_properties.border_color.type == TC_LINESTYLE) {
747 struct lp_style_type style;
748 int side = pm3d_side( &pointsA[i], &pointsA[i1], &pointsB[ii]);
749 lp_use_properties(&style, side < 0 ? cb1 + 1 : cb1);
750 cb1 = cb2 = cb3 = cb4 = style.pm3d_color.lt;
751 }
752 } else {
753 cb1 = pointsA[i].z;
754 cb2 = pointsA[i1].z;
755 cb3 = pointsB[ii].z;
756 cb4 = pointsB[ii1].z;
757 }
758 /* Fancy averages of RGB color make no sense */
759 if (color_from_rgbvar) {
760 unsigned int r, g, b, a;
761 unsigned int u1 = cb1;
762 unsigned int u2 = cb2;
763 unsigned int u3 = cb3;
764 unsigned int u4 = cb4;
765 switch (pm3d.which_corner_color) {
766 default:
767 r = (u1&0xff0000) + (u2&0xff0000) + (u3&0xff0000) + (u4&0xff0000);
768 g = (u1&0xff00) + (u2&0xff00) + (u3&0xff00) + (u4&0xff00);
769 b = (u1&0xff) + (u2&0xff) + (u3&0xff) + (u4&0xff);
770 avgC = ((r>>2)&0xff0000) + ((g>>2)&0xff00) + ((b>>2)&0xff);
771 a = ((u1>>24)&0xff) + ((u2>>24)&0xff) + ((u3>>24)&0xff) + ((u4>>24)&0xff);
772 avgC += (a<<22)&0xff000000;
773 break;
774 case PM3D_WHICHCORNER_C1: avgC = cb1; break;
775 case PM3D_WHICHCORNER_C2: avgC = cb2; break;
776 case PM3D_WHICHCORNER_C3: avgC = cb3; break;
777 case PM3D_WHICHCORNER_C4: avgC = cb4; break;
778 }
779
780 /* But many different averages are possible for gray values */
781 } else {
782 switch (pm3d.which_corner_color) {
783 default:
784 case PM3D_WHICHCORNER_MEAN: avgC = (cb1 + cb2 + cb3 + cb4) * 0.25; break;
785 case PM3D_WHICHCORNER_GEOMEAN: avgC = geomean4(cb1, cb2, cb3, cb4); break;
786 case PM3D_WHICHCORNER_HARMEAN: avgC = harmean4(cb1, cb2, cb3, cb4); break;
787 case PM3D_WHICHCORNER_MEDIAN: avgC = median4(cb1, cb2, cb3, cb4); break;
788 case PM3D_WHICHCORNER_MIN: avgC = minimum4(cb1, cb2, cb3, cb4); break;
789 case PM3D_WHICHCORNER_MAX: avgC = maximum4(cb1, cb2, cb3, cb4); break;
790 case PM3D_WHICHCORNER_RMS: avgC = rms4(cb1, cb2, cb3, cb4); break;
791 case PM3D_WHICHCORNER_C1: avgC = cb1; break;
792 case PM3D_WHICHCORNER_C2: avgC = cb2; break;
793 case PM3D_WHICHCORNER_C3: avgC = cb3; break;
794 case PM3D_WHICHCORNER_C4: avgC = cb4; break;
795 }
796 }
797
798 /* The value is out of range, but we didn't figure it out until now */
799 if (isnan(avgC))
800 continue;
801
802 /* Option to not drawn quadrangles with cb out of range */
803 if (pm3d.no_clipcb && (avgC > CB_AXIS.max || avgC < CB_AXIS.min))
804 continue;
805
806 if (color_from_rgbvar) /* we were given an RGB color */
807 gray = avgC;
808 else /* transform z value to gray, i.e. to interval [0,1] */
809 gray = cb2gray(avgC);
810
811 /* apply lighting model */
812 if (pm3d_shade.strength > 0) {
813 if (at_which_z == PM3D_AT_SURFACE)
814 gray = apply_lighting_model( &pointsA[i], &pointsA[i1],
815 &pointsB[ii], &pointsB[ii1], gray);
816 /* Don't apply lighting model to TOP/BOTTOM projections */
817 /* but convert from floating point 0<gray<1 to RGB color */
818 /* since that is what would have been returned from the */
819 /* lighting code. */
820 else if (!color_from_rgbvar) {
821 rgb255_color temp;
822 rgb255maxcolors_from_gray(gray, &temp);
823 gray = (long)((temp.r << 16) + (temp.g << 8) + (temp.b));
824 }
825 }
826
827 /* set the color */
828 if (pm3d.direction != PM3D_DEPTH) {
829 if (color_from_rgbvar || pm3d_shade.strength > 0)
830 set_rgbcolor_var((unsigned int)gray);
831 else
832 set_color(gray);
833 }
834 }
835 }
836
837 corners[0].x = pointsA[i].x;
838 corners[0].y = pointsA[i].y;
839 corners[1].x = pointsB[ii].x;
840 corners[1].y = pointsB[ii].y;
841 corners[2].x = pointsB[ii1].x;
842 corners[2].y = pointsB[ii1].y;
843 corners[3].x = pointsA[i1].x;
844 corners[3].y = pointsA[i1].y;
845
846 if (interp_i > 1 || interp_j > 1 || at_which_z == PM3D_AT_SURFACE) {
847 corners[0].z = pointsA[i].z;
848 corners[1].z = pointsB[ii].z;
849 corners[2].z = pointsB[ii1].z;
850 corners[3].z = pointsA[i1].z;
851 if (color_from_column) {
852 corners[0].c = pointsA[i].CRD_COLOR;
853 corners[1].c = pointsB[ii].CRD_COLOR;
854 corners[2].c = pointsB[ii1].CRD_COLOR;
855 corners[3].c = pointsA[i1].CRD_COLOR;
856 }
857 }
858 if (interp_i > 1 || interp_j > 1) {
859 /* Interpolation is enabled.
860 * interp_i is the # of points along scan lines
861 * interp_j is the # of points between scan lines
862 * Algorithm is to first sample i points along the scan lines
863 * defined by corners[3],corners[0] and corners[2],corners[1]. */
864 int j1;
865 for (i1 = 0; i1 <= interp_i; i1++) {
866 bl_point[i1][0].x =
867 ((corners[3].x - corners[0].x) / interp_i) * i1 + corners[0].x;
868 bl_point[i1][interp_j].x =
869 ((corners[2].x - corners[1].x) / interp_i) * i1 + corners[1].x;
870 bl_point[i1][0].y =
871 ((corners[3].y - corners[0].y) / interp_i) * i1 + corners[0].y;
872 bl_point[i1][interp_j].y =
873 ((corners[2].y - corners[1].y) / interp_i) * i1 + corners[1].y;
874 bl_point[i1][0].z =
875 ((corners[3].z - corners[0].z) / interp_i) * i1 + corners[0].z;
876 bl_point[i1][interp_j].z =
877 ((corners[2].z - corners[1].z) / interp_i) * i1 + corners[1].z;
878 if (color_from_column) {
879 bl_point[i1][0].c =
880 ((corners[3].c - corners[0].c) / interp_i) * i1 + corners[0].c;
881 bl_point[i1][interp_j].c =
882 ((corners[2].c - corners[1].c) / interp_i) * i1 + corners[1].c;
883 }
884 /* Next we sample j points between each of the new points
885 * created in the previous step (this samples between
886 * scan lines) in the same manner. */
887 for (j1 = 1; j1 < interp_j; j1++) {
888 bl_point[i1][j1].x =
889 ((bl_point[i1][interp_j].x - bl_point[i1][0].x) / interp_j) *
890 j1 + bl_point[i1][0].x;
891 bl_point[i1][j1].y =
892 ((bl_point[i1][interp_j].y - bl_point[i1][0].y) / interp_j) *
893 j1 + bl_point[i1][0].y;
894 bl_point[i1][j1].z =
895 ((bl_point[i1][interp_j].z - bl_point[i1][0].z) / interp_j) *
896 j1 + bl_point[i1][0].z;
897 if (color_from_column)
898 bl_point[i1][j1].c =
899 ((bl_point[i1][interp_j].c - bl_point[i1][0].c) / interp_j) *
900 j1 + bl_point[i1][0].c;
901 }
902 }
903 /* Once all points are created, move them into an appropriate
904 * structure and call set_color on each to retrieve the
905 * correct color mapping for this new sub-sampled quadrangle. */
906 for (i1 = 0; i1 < interp_i; i1++) {
907 for (j1 = 0; j1 < interp_j; j1++) {
908 corners[0].x = bl_point[i1][j1].x;
909 corners[0].y = bl_point[i1][j1].y;
910 corners[0].z = bl_point[i1][j1].z;
911 corners[1].x = bl_point[i1+1][j1].x;
912 corners[1].y = bl_point[i1+1][j1].y;
913 corners[1].z = bl_point[i1+1][j1].z;
914 corners[2].x = bl_point[i1+1][j1+1].x;
915 corners[2].y = bl_point[i1+1][j1+1].y;
916 corners[2].z = bl_point[i1+1][j1+1].z;
917 corners[3].x = bl_point[i1][j1+1].x;
918 corners[3].y = bl_point[i1][j1+1].y;
919 corners[3].z = bl_point[i1][j1+1].z;
920 if (color_from_column) {
921 corners[0].c = bl_point[i1][j1].c;
922 corners[1].c = bl_point[i1+1][j1].c;
923 corners[2].c = bl_point[i1+1][j1+1].c;
924 corners[3].c = bl_point[i1][j1+1].c;
925 }
926
927 FPRINTF((stderr,"(%g,%g),(%g,%g),(%g,%g),(%g,%g)\n",
928 corners[0].x, corners[0].y, corners[1].x, corners[1].y,
929 corners[2].x, corners[2].y, corners[3].x, corners[3].y));
930
931 /* If the colors are given separately, we already loaded them above */
932 if (color_from_column) {
933 cb1 = corners[0].c;
934 cb2 = corners[1].c;
935 cb3 = corners[2].c;
936 cb4 = corners[3].c;
937 } else {
938 cb1 = corners[0].z;
939 cb2 = corners[1].z;
940 cb3 = corners[2].z;
941 cb4 = corners[3].z;
942 }
943 switch (pm3d.which_corner_color) {
944 default:
945 case PM3D_WHICHCORNER_MEAN: avgC = (cb1 + cb2 + cb3 + cb4) * 0.25; break;
946 case PM3D_WHICHCORNER_GEOMEAN: avgC = geomean4(cb1, cb2, cb3, cb4); break;
947 case PM3D_WHICHCORNER_HARMEAN: avgC = harmean4(cb1, cb2, cb3, cb4); break;
948 case PM3D_WHICHCORNER_MEDIAN: avgC = median4(cb1, cb2, cb3, cb4); break;
949 case PM3D_WHICHCORNER_MIN: avgC = minimum4(cb1, cb2, cb3, cb4); break;
950 case PM3D_WHICHCORNER_MAX: avgC = maximum4(cb1, cb2, cb3, cb4); break;
951 case PM3D_WHICHCORNER_RMS: avgC = rms4(cb1, cb2, cb3, cb4); break;
952 case PM3D_WHICHCORNER_C1: avgC = cb1; break;
953 case PM3D_WHICHCORNER_C2: avgC = cb2; break;
954 case PM3D_WHICHCORNER_C3: avgC = cb3; break;
955 case PM3D_WHICHCORNER_C4: avgC = cb4; break;
956 }
957
958 if (color_from_fillcolor) {
959 /* color is set by "fc <rgbval>" */
960 gray = this_plot->fill_properties.border_color.lt;
961 /* EXPERIMENTAL
962 * pm3d fc linestyle N generates
963 * top/bottom color difference as with hidden3d
964 */
965 if (this_plot->fill_properties.border_color.type == TC_LINESTYLE) {
966 struct lp_style_type style;
967 int side = pm3d_side( &pointsA[i], &pointsB[ii], &pointsB[ii1]);
968 lp_use_properties(&style, side < 0 ? gray + 1 : gray);
969 gray = style.pm3d_color.lt;
970 }
971 } else if (color_from_rgbvar) {
972 /* we were given an explicit color */
973 gray = avgC;
974 } else {
975 /* clip if out of range */
976 if (pm3d.no_clipcb && (avgC < CB_AXIS.min || avgC > CB_AXIS.max))
977 continue;
978 /* transform z value to gray, i.e. to interval [0,1] */
979 gray = cb2gray(avgC);
980 }
981
982 /* apply lighting model */
983 if (pm3d_shade.strength > 0) {
984 /* FIXME: coordinate->quadrangle->coordinate seems crazy */
985 coordinate corcorners[4];
986 int i;
987 for (i=0; i<4; i++) {
988 corcorners[i].x = corners[i].x;
989 corcorners[i].y = corners[i].y;
990 corcorners[i].z = corners[i].z;
991 }
992
993 if (at_which_z == PM3D_AT_SURFACE)
994 gray = apply_lighting_model( &corcorners[0], &corcorners[1],
995 &corcorners[2], &corcorners[3], gray);
996 /* Don't apply lighting model to TOP/BOTTOM projections */
997 /* but convert from floating point 0<gray<1 to RGB color */
998 /* since that is what would have been returned from the */
999 /* lighting code. */
1000 else if (!color_from_rgbvar) {
1001 rgb255_color temp;
1002 rgb255maxcolors_from_gray(gray, &temp);
1003 gray = (long)((temp.r << 16) + (temp.g << 8) + (temp.b));
1004 }
1005 }
1006
1007 if (pm3d.direction == PM3D_DEPTH) {
1008 /* copy quadrangle */
1009 quadrangle* qp = quadrangles + current_quadrangle;
1010 memcpy(qp->vertex.corners, corners, 4 * sizeof (gpdPoint));
1011 if (color_from_rgbvar) {
1012 qp->gray = PM3D_USE_RGB_COLOR_INSTEAD_OF_GRAY;
1013 qp->qcolor.rgb_color = (unsigned int)gray;
1014 } else {
1015 qp->gray = gray;
1016 qp->qcolor.colorspec = &this_plot->lp_properties.pm3d_color;
1017 }
1018 qp->fillstyle = plot_fillstyle;
1019 qp->type = QUAD_TYPE_NORMAL;
1020 current_quadrangle++;
1021 } else {
1022 if (pm3d_shade.strength > 0 || color_from_rgbvar)
1023 set_rgbcolor_var((unsigned int)gray);
1024 else
1025 set_color(gray);
1026 if (at_which_z == PM3D_AT_BASE)
1027 corners[0].z = corners[1].z = corners[2].z = corners[3].z = base_z;
1028 else if (at_which_z == PM3D_AT_TOP)
1029 corners[0].z = corners[1].z = corners[2].z = corners[3].z = ceiling_z;
1030 filled_polygon(corners, plot_fillstyle, 4);
1031 }
1032 }
1033 }
1034 } else { /* thus (interp_i == 1 && interp_j == 1) */
1035 if (pm3d.direction != PM3D_DEPTH) {
1036 filled_polygon(corners, plot_fillstyle, 4);
1037 } else {
1038 /* copy quadrangle */
1039 quadrangle* qp = quadrangles + current_quadrangle;
1040 memcpy(qp->vertex.corners, corners, 4 * sizeof (gpdPoint));
1041 if (color_from_rgbvar) {
1042 qp->gray = PM3D_USE_RGB_COLOR_INSTEAD_OF_GRAY;
1043 qp->qcolor.rgb_color = (unsigned int)gray;
1044 } else {
1045 qp->gray = gray;
1046 qp->qcolor.colorspec = &this_plot->lp_properties.pm3d_color;
1047 }
1048 qp->fillstyle = plot_fillstyle;
1049 qp->type = QUAD_TYPE_NORMAL;
1050 current_quadrangle++;
1051 }
1052 } /* interpolate between points */
1053 } /* loop quadrangles over points of two subsequent scans */
1054 } /* loop over scans */
1055
1056 if (bl_point) {
1057 for (i1 = 0; i1 <= interp_i; i1++)
1058 free(bl_point[i1]);
1059 free(bl_point);
1060 }
1061 /* free memory allocated by scan_array */
1062 free(scan_array);
1063
1064 /* reset local flags */
1065 pm3d_plot_at = 0;
1066
1067 /* for pm3dCompress.awk and pm3dConvertToImage.awk */
1068 if (pm3d.direction != PM3D_DEPTH)
1069 term->layer(TERM_LAYER_END_PM3D_MAP);
1070 } /* end of pm3d splotting mode */
1071
1072
1073 /*
1074 * unset pm3d for the reset command
1075 */
1076 void
pm3d_reset()1077 pm3d_reset()
1078 {
1079 strcpy(pm3d.where, "s");
1080 pm3d_plot_at = 0;
1081 pm3d.flush = PM3D_FLUSH_BEGIN;
1082 pm3d.ftriangles = 0;
1083 pm3d.clip = PM3D_CLIP_Z; /* prior to Dec 2019 default was PM3D_CLIP_4IN */
1084 pm3d.no_clipcb = FALSE;
1085 pm3d.direction = PM3D_SCANS_AUTOMATIC;
1086 pm3d.base_sort = FALSE;
1087 pm3d.zmean_sort = TRUE; /* DEBUG: prior to Oct 2019 default sort used zmax */
1088 pm3d.implicit = PM3D_EXPLICIT;
1089 pm3d.which_corner_color = PM3D_WHICHCORNER_MEAN;
1090 pm3d.interp_i = 1;
1091 pm3d.interp_j = 1;
1092 pm3d.border = default_pm3d_border;
1093 pm3d.border.l_type = LT_NODRAW;
1094
1095 pm3d_shade.strength = 0.0;
1096 pm3d_shade.spec = 0.0;
1097 pm3d_shade.fixed = TRUE;
1098 }
1099
1100
1101 /*
1102 * Draw (one) PM3D color surface.
1103 */
1104 void
pm3d_draw_one(struct surface_points * plot)1105 pm3d_draw_one(struct surface_points *plot)
1106 {
1107 int i = 0;
1108 char *where = plot->pm3d_where[0] ? plot->pm3d_where : pm3d.where;
1109 /* Draw either at 'where' option of the given surface or at pm3d.where
1110 * global option. */
1111
1112 if (!where[0])
1113 return;
1114
1115 /* Initialize lighting model */
1116 if (pm3d_shade.strength > 0)
1117 pm3d_init_lighting_model();
1118
1119 for (; where[i]; i++) {
1120 pm3d_plot(plot, where[i]);
1121 }
1122 }
1123
1124
1125 /*
1126 * Add one pm3d quadrangle to the mix.
1127 * Called by zerrorfill() and by plot3d_boxes().
1128 * Also called by vplot_isosurface().
1129 */
1130 void
pm3d_add_quadrangle(struct surface_points * plot,gpdPoint corners[4])1131 pm3d_add_quadrangle(struct surface_points *plot, gpdPoint corners[4])
1132 {
1133 pm3d_add_polygon(plot, corners, 4);
1134 }
1135
1136 /*
1137 * The general case.
1138 * (plot == NULL) if we were called from do_polygon().
1139 */
1140 void
pm3d_add_polygon(struct surface_points * plot,gpdPoint corners[4],int vertices)1141 pm3d_add_polygon(struct surface_points *plot, gpdPoint corners[4], int vertices)
1142 {
1143 quadrangle *q;
1144
1145 /* FIXME: I have no idea how to estimate the number of facets for an isosurface */
1146 if (!plot || (plot->plot_style == ISOSURFACE)) {
1147 if (allocated_quadrangles < current_quadrangle + 100) {
1148 allocated_quadrangles += 1000.;
1149 quadrangles = gp_realloc(quadrangles,
1150 allocated_quadrangles * sizeof(quadrangle), "pm3d_add_quadrangle");
1151 }
1152 } else
1153
1154 if (allocated_quadrangles < current_quadrangle + plot->iso_crvs->p_count) {
1155 allocated_quadrangles += 2 * plot->iso_crvs->p_count;
1156 quadrangles = gp_realloc(quadrangles,
1157 allocated_quadrangles * sizeof(quadrangle), "pm3d_add_quadrangle");
1158 }
1159 q = quadrangles + current_quadrangle++;
1160 memcpy(q->vertex.corners, corners, 4*sizeof(gpdPoint));
1161 if (plot)
1162 q->fillstyle = style_from_fill(&plot->fill_properties);
1163 else
1164 q->fillstyle = 0;
1165
1166 /* For triangles and normal quadrangles, the vertices are stored in
1167 * q->vertex.corners. For larger polygons we store them in external array
1168 * polygonlist and keep only an index into the array q->vertex.array.
1169 */
1170 q->type = QUAD_TYPE_NORMAL;
1171 if (corners[3].x == corners[2].x && corners[3].y == corners[2].y
1172 && corners[3].z == corners[2].z)
1173 q->type = QUAD_TYPE_TRIANGLE;
1174 if (vertices > 4) {
1175 gpdPoint *save_corners = get_polygon(vertices);
1176 q->vertex.array_index = current_polygon;
1177 q->type = QUAD_TYPE_LARGEPOLYGON;
1178 memcpy(save_corners, corners, vertices * sizeof(gpdPoint));
1179 save_corners[2].c = vertices;
1180 }
1181
1182 if (!plot) {
1183 /* This quadrangle came from "set object polygon" rather than "splot with pm3d" */
1184 if (corners[0].c == LT_BACKGROUND) {
1185 q->gray = PM3D_USE_BACKGROUND_INSTEAD_OF_GRAY;
1186 } else {
1187 q->qcolor.rgb_color = corners[0].c;
1188 q->gray = PM3D_USE_RGB_COLOR_INSTEAD_OF_GRAY;
1189 }
1190 q->fillstyle = corners[1].c;
1191
1192 } else if (plot->pm3d_color_from_column) {
1193 /* FIXME: color_from_rgbvar need only be set once per plot */
1194 /* This is the usual path for 'splot with boxes' */
1195 color_from_rgbvar = TRUE;
1196 if (pm3d_shade.strength > 0) {
1197 q->gray = plot->lp_properties.pm3d_color.lt;
1198 illuminate_one_quadrangle(q);
1199 } else {
1200 q->qcolor.rgb_color = plot->lp_properties.pm3d_color.lt;
1201 q->gray = PM3D_USE_RGB_COLOR_INSTEAD_OF_GRAY;
1202 }
1203
1204 } else if (plot->lp_properties.pm3d_color.type == TC_Z) {
1205 /* This is a special case for 'splot with boxes lc palette z' */
1206 q->gray = cb2gray(corners[1].z);
1207 color_from_rgbvar = FALSE;
1208 if (pm3d_shade.strength > 0)
1209 illuminate_one_quadrangle(q);
1210
1211 } else if (plot->plot_style == ISOSURFACE
1212 || plot->plot_style == POLYGONS) {
1213
1214 int rgb_color = corners[0].c;
1215 if (corners[0].c == LT_BACKGROUND)
1216 q->gray = PM3D_USE_BACKGROUND_INSTEAD_OF_GRAY;
1217 else
1218 q->gray = PM3D_USE_RGB_COLOR_INSTEAD_OF_GRAY;
1219
1220 if (plot->plot_style == ISOSURFACE && isosurface_options.inside_offset > 0) {
1221 struct lp_style_type style;
1222 struct coordinate v[3];
1223 int i;
1224 for (i=0; i<3; i++) {
1225 v[i].x = corners[i].x;
1226 v[i].y = corners[i].y;
1227 v[i].z = corners[i].z;
1228 }
1229 i = plot->hidden3d_top_linetype + 1;
1230 if (pm3d_side( &v[0], &v[1], &v[2] ) < 0)
1231 i += isosurface_options.inside_offset;
1232 lp_use_properties(&style, i);
1233 rgb_color = style.pm3d_color.lt;
1234 }
1235 q->qcolor.rgb_color = rgb_color;
1236 if (pm3d_shade.strength > 0) {
1237 q->gray = rgb_color;
1238 color_from_rgbvar = TRUE;
1239 illuminate_one_quadrangle(q);
1240 }
1241
1242 } else {
1243 /* This is the usual [only?] path for 'splot with zerror' */
1244 q->qcolor.colorspec = &plot->fill_properties.border_color;
1245 q->gray = PM3D_USE_COLORSPEC_INSTEAD_OF_GRAY;
1246 }
1247 }
1248
1249
1250
1251 /* Display an error message for the routine get_pm3d_at_option() below.
1252 */
1253 static void
pm3d_option_at_error()1254 pm3d_option_at_error()
1255 {
1256 int_error(c_token, "\
1257 parameter to `pm3d at` requires combination of up to 6 characters b,s,t\n\
1258 \t(drawing at bottom, surface, top)");
1259 }
1260
1261
1262 /* Read the option for 'pm3d at' command.
1263 * Used by 'set pm3d at ...' or by 'splot ... with pm3d at ...'.
1264 * If no option given, then returns empty string, otherwise copied there.
1265 * The string is unchanged on error, and 1 is returned.
1266 * On success, 0 is returned.
1267 */
1268 int
get_pm3d_at_option(char * pm3d_where)1269 get_pm3d_at_option(char *pm3d_where)
1270 {
1271 char* c;
1272
1273 if (END_OF_COMMAND || token[c_token].length >= sizeof(pm3d.where)) {
1274 pm3d_option_at_error();
1275 return 1;
1276 }
1277 memcpy(pm3d_where, gp_input_line + token[c_token].start_index,
1278 token[c_token].length);
1279 pm3d_where[ token[c_token].length ] = 0;
1280 /* verify the parameter */
1281 for (c = pm3d_where; *c; c++) {
1282 if (*c != PM3D_AT_BASE && *c != PM3D_AT_TOP && *c != PM3D_AT_SURFACE) {
1283 pm3d_option_at_error();
1284 return 1;
1285 }
1286 }
1287 c_token++;
1288 return 0;
1289 }
1290
1291 /* Set flag plot_has_palette to TRUE if there is any element on the graph
1292 * which requires palette of continuous colors.
1293 */
1294 void
set_plot_with_palette(int plot_num,int plot_mode)1295 set_plot_with_palette(int plot_num, int plot_mode)
1296 {
1297 struct surface_points *this_3dplot = first_3dplot;
1298 struct curve_points *this_2dplot = first_plot;
1299 int surface = 0;
1300 struct text_label *this_label = first_label;
1301 struct object *this_object;
1302
1303 plot_has_palette = TRUE;
1304 /* Is pm3d switched on globally? */
1305 if (pm3d.implicit == PM3D_IMPLICIT)
1306 return;
1307
1308 /* Check 2D plots */
1309 if (plot_mode == MODE_PLOT) {
1310 while (this_2dplot) {
1311 if (this_2dplot->plot_style == IMAGE)
1312 return;
1313 if (this_2dplot->lp_properties.pm3d_color.type == TC_CB
1314 || this_2dplot->lp_properties.pm3d_color.type == TC_FRAC
1315 || this_2dplot->lp_properties.pm3d_color.type == TC_Z)
1316 return;
1317 if (this_2dplot->labels
1318 && (this_2dplot->labels->textcolor.type == TC_CB
1319 || this_2dplot->labels->textcolor.type == TC_FRAC
1320 || this_2dplot->labels->textcolor.type == TC_Z))
1321 return;
1322 this_2dplot = this_2dplot->next;
1323 }
1324 }
1325
1326 /* Check 3D plots */
1327 if (plot_mode == MODE_SPLOT) {
1328 /* Any surface 'with pm3d', 'with image' or 'with line|dot palette'? */
1329 while (surface++ < plot_num) {
1330 int type;
1331 if (this_3dplot->plot_style == PM3DSURFACE)
1332 return;
1333 if (this_3dplot->plot_style == IMAGE)
1334 return;
1335
1336 type = this_3dplot->lp_properties.pm3d_color.type;
1337 if (type == TC_LT || type == TC_LINESTYLE || type == TC_RGB)
1338 ; /* don't return yet */
1339 else
1340 /* TC_DEFAULT: splot x with line|lp|dot palette */
1341 return;
1342
1343 if (this_3dplot->labels &&
1344 this_3dplot->labels->textcolor.type >= TC_CB)
1345 return;
1346 this_3dplot = this_3dplot->next_sp;
1347 }
1348 }
1349
1350 #define TC_USES_PALETTE(tctype) (tctype==TC_Z) || (tctype==TC_CB) || (tctype==TC_FRAC)
1351
1352 /* Any label with 'textcolor palette'? */
1353 for (; this_label != NULL; this_label = this_label->next) {
1354 if (TC_USES_PALETTE(this_label->textcolor.type))
1355 return;
1356 }
1357 /* Any of title, xlabel, ylabel, zlabel, ... with 'textcolor palette'? */
1358 if (TC_USES_PALETTE(title.textcolor.type)) return;
1359 if (TC_USES_PALETTE(axis_array[FIRST_X_AXIS].label.textcolor.type)) return;
1360 if (TC_USES_PALETTE(axis_array[FIRST_Y_AXIS].label.textcolor.type)) return;
1361 if (TC_USES_PALETTE(axis_array[SECOND_X_AXIS].label.textcolor.type)) return;
1362 if (TC_USES_PALETTE(axis_array[SECOND_Y_AXIS].label.textcolor.type)) return;
1363 if (plot_mode == MODE_SPLOT)
1364 if (TC_USES_PALETTE(axis_array[FIRST_Z_AXIS].label.textcolor.type)) return;
1365 if (TC_USES_PALETTE(axis_array[COLOR_AXIS].label.textcolor.type)) return;
1366 for (this_object = first_object; this_object != NULL; this_object = this_object->next) {
1367 if (TC_USES_PALETTE(this_object->lp_properties.pm3d_color.type))
1368 return;
1369 }
1370
1371 #undef TC_USES_PALETTE
1372
1373 /* Palette with continuous colors is not used. */
1374 plot_has_palette = FALSE; /* otherwise it stays TRUE */
1375 }
1376
1377 TBOOLEAN
is_plot_with_palette()1378 is_plot_with_palette()
1379 {
1380 return plot_has_palette;
1381 }
1382
1383 TBOOLEAN
is_plot_with_colorbox()1384 is_plot_with_colorbox()
1385 {
1386 return plot_has_palette && (color_box.where != SMCOLOR_BOX_NO);
1387 }
1388
1389 /*
1390 * Must be called before trying to apply lighting model
1391 */
1392 void
pm3d_init_lighting_model()1393 pm3d_init_lighting_model()
1394 {
1395 light[0] = cos(-DEG2RAD*pm3d_shade.rot_x)*cos(-(DEG2RAD*pm3d_shade.rot_z+90));
1396 light[2] = cos(-DEG2RAD*pm3d_shade.rot_x)*sin(-(DEG2RAD*pm3d_shade.rot_z+90));
1397 light[1] = sin(-DEG2RAD*pm3d_shade.rot_x);
1398 }
1399
1400 /*
1401 * Layer on layer of coordinate conventions
1402 */
1403 static void
illuminate_one_quadrangle(quadrangle * q)1404 illuminate_one_quadrangle( quadrangle *q )
1405 {
1406 struct coordinate c1, c2, c3, c4;
1407 vertex vtmp;
1408
1409 map3d_xyz(q->vertex.corners[0].x, q->vertex.corners[0].y, q->vertex.corners[0].z, &vtmp);
1410 c1.x = vtmp.x; c1.y = vtmp.y; c1.z = vtmp.z;
1411 map3d_xyz(q->vertex.corners[1].x, q->vertex.corners[1].y, q->vertex.corners[1].z, &vtmp);
1412 c2.x = vtmp.x; c2.y = vtmp.y; c2.z = vtmp.z;
1413 map3d_xyz(q->vertex.corners[2].x, q->vertex.corners[2].y, q->vertex.corners[2].z, &vtmp);
1414 c3.x = vtmp.x; c3.y = vtmp.y; c3.z = vtmp.z;
1415 map3d_xyz(q->vertex.corners[3].x, q->vertex.corners[3].y, q->vertex.corners[3].z, &vtmp);
1416 c4.x = vtmp.x; c4.y = vtmp.y; c4.z = vtmp.z;
1417 q->gray = apply_lighting_model( &c1, &c2, &c3, &c4, q->gray);
1418 }
1419
1420 /*
1421 * Adjust current RGB color based on pm3d lighting model.
1422 * Jan 2019: preserve alpha channel
1423 * This isn't quite right because specular highlights should
1424 * not be affected by transparency.
1425 */
1426 int
apply_lighting_model(struct coordinate * v0,struct coordinate * v1,struct coordinate * v2,struct coordinate * v3,double gray)1427 apply_lighting_model( struct coordinate *v0, struct coordinate *v1,
1428 struct coordinate *v2, struct coordinate *v3,
1429 double gray )
1430 {
1431 double normal[3];
1432 double normal1[3];
1433 double reflect[3];
1434 double t;
1435 double phi;
1436 double psi;
1437 unsigned int rgb;
1438 unsigned int alpha = 0;
1439 rgb_color color;
1440 double r, g, b, tmp_r, tmp_g, tmp_b;
1441 double dot_prod, shade_fact, spec_fact;
1442
1443 if (color_from_rgbvar) {
1444 rgb = gray;
1445 r = (double)((rgb >> 16) & 0xFF) / 255.;
1446 g = (double)((rgb >> 8) & 0xFF) / 255.;
1447 b = (double)((rgb ) & 0xFF) / 255.;
1448 alpha = rgb & 0xff000000;
1449 } else {
1450 rgb1_from_gray(gray, &color);
1451 r = color.r;
1452 g = color.g;
1453 b = color.b;
1454 }
1455
1456 psi = -DEG2RAD*(surface_rot_z);
1457 phi = -DEG2RAD*(surface_rot_x);
1458
1459 normal[0] = (v1->y-v0->y)*(v2->z-v0->z)*yscale3d*zscale3d
1460 - (v1->z-v0->z)*(v2->y-v0->y)*yscale3d*zscale3d;
1461 normal[1] = (v1->z-v0->z)*(v2->x-v0->x)*xscale3d*zscale3d
1462 - (v1->x-v0->x)*(v2->z-v0->z)*xscale3d*zscale3d;
1463 normal[2] = (v1->x-v0->x)*(v2->y-v0->y)*xscale3d*yscale3d
1464 - (v1->y-v0->y)*(v2->x-v0->x)*xscale3d*yscale3d ;
1465
1466 t = sqrt( normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2] );
1467
1468 /* Trap and handle degenerate case of two identical vertices.
1469 * Aug 2020
1470 * FIXME: The problem case that revealed this problem always involved
1471 * v2 == (v0 or v1) but some unknown example might instead yield
1472 * v0 == v1, which is not addressed here.
1473 */
1474 if (t < 1.e-12) {
1475 if (v2 == v3) /* 2nd try; give up and return original color */
1476 return (color_from_rgbvar)
1477 ? gray
1478 : ((unsigned char)(r*255.) << 16)
1479 + ((unsigned char)(g*255.) << 8)
1480 + ((unsigned char)(b*255.));
1481 else
1482 return apply_lighting_model(v0, v1, v3, v3, gray);
1483 }
1484
1485 normal[0] /= t;
1486 normal[1] /= t;
1487 normal[2] /= t;
1488
1489 /* Correct for the view angle so that the illumination is "fixed" with */
1490 /* respect to the viewer rather than rotating with the surface. */
1491 if (pm3d_shade.fixed) {
1492 normal1[0] = cos(psi)*normal[0] - sin(psi)*normal[1] + 0*normal[2];
1493 normal1[1] = sin(psi)*normal[0] + cos(psi)*normal[1] + 0*normal[2];
1494 normal1[2] = 0*normal[0] + 0*normal[1] + 1*normal[2];
1495
1496 normal[0] = 1*normal1[0] + 0*normal1[1] + 0*normal1[2];
1497 normal[1] = 0*normal1[0] + cos(phi)*normal1[1] - sin(phi)*normal1[2];
1498 normal[2] = 0*normal1[0] + sin(phi)*normal1[1] + cos(phi)*normal1[2];
1499 }
1500
1501 if (normal[2] < 0.0) {
1502 normal[0] *= -1.0;
1503 normal[1] *= -1.0;
1504 normal[2] *= -1.0;
1505 }
1506
1507 dot_prod = normal[0]*light[0] + normal[1]*light[1] + normal[2]*light[2];
1508 shade_fact = (dot_prod < 0) ? -dot_prod : 0;
1509
1510 tmp_r = r*(pm3d_shade.ambient-pm3d_shade.strength+shade_fact*pm3d_shade.strength);
1511 tmp_g = g*(pm3d_shade.ambient-pm3d_shade.strength+shade_fact*pm3d_shade.strength);
1512 tmp_b = b*(pm3d_shade.ambient-pm3d_shade.strength+shade_fact*pm3d_shade.strength);
1513
1514 /* Specular highlighting */
1515 if (pm3d_shade.spec > 0.0) {
1516
1517 reflect[0] = -light[0]+2*dot_prod*normal[0];
1518 reflect[1] = -light[1]+2*dot_prod*normal[1];
1519 reflect[2] = -light[2]+2*dot_prod*normal[2];
1520 t = sqrt( reflect[0]*reflect[0] + reflect[1]*reflect[1] + reflect[2]*reflect[2] );
1521 reflect[0] /= t;
1522 reflect[1] /= t;
1523 reflect[2] /= t;
1524
1525 dot_prod = -reflect[2];
1526
1527 /* old-style Phong equation; no need for bells or whistles */
1528 spec_fact = pow(fabs(dot_prod), pm3d_shade.Phong);
1529
1530 /* White spotlight from direction phi,psi */
1531 if (dot_prod > 0) {
1532 tmp_r += pm3d_shade.spec * spec_fact;
1533 tmp_g += pm3d_shade.spec * spec_fact;
1534 tmp_b += pm3d_shade.spec * spec_fact;
1535 }
1536 /* EXPERIMENTAL: Red spotlight from underneath */
1537 if (dot_prod < 0 && pm3d_shade.spec2 > 0) {
1538 tmp_r += pm3d_shade.spec2 * spec_fact;
1539 }
1540 }
1541
1542 tmp_r = clip_to_01(tmp_r);
1543 tmp_g = clip_to_01(tmp_g);
1544 tmp_b = clip_to_01(tmp_b);
1545
1546 rgb = ((unsigned char)((tmp_r)*255.) << 16)
1547 + ((unsigned char)((tmp_g)*255.) << 8)
1548 + ((unsigned char)((tmp_b)*255.));
1549
1550 /* restore alpha value if there was one */
1551 rgb |= alpha;
1552
1553 return rgb;
1554 }
1555
1556
1557 /* The pm3d code works with gpdPoint data structures (double: x,y,z,color)
1558 * term->filled_polygon want gpiPoint data (int: x,y,style).
1559 * This routine converts from gpdPoint to gpiPoint
1560 */
1561 static void
filled_polygon(gpdPoint * corners,int fillstyle,int nv)1562 filled_polygon(gpdPoint *corners, int fillstyle, int nv)
1563 {
1564 int i;
1565 double x, y;
1566
1567 /* For normal pm3d surfaces and tessellation the constituent polygons
1568 * have a small number of vertices, usually 4.
1569 * However generalized polygons like cartographic areas could have a
1570 * vastly larger number of vertices, so we allow for dynamic reallocation.
1571 */
1572 static int max_vertices = 0;
1573 static gpiPoint *icorners = NULL;
1574 static gpiPoint *ocorners = NULL;
1575 static gpdPoint *clipcorners = NULL;
1576 if (nv > max_vertices) {
1577 max_vertices = nv;
1578 icorners = gp_realloc( icorners, (2*max_vertices) * sizeof(gpiPoint), "filled_polygon");
1579 ocorners = gp_realloc( ocorners, (2*max_vertices) * sizeof(gpiPoint), "filled_polygon");
1580 clipcorners = gp_realloc( clipcorners, (2*max_vertices) * sizeof(gpdPoint), "filled_polygon");
1581 }
1582
1583 if ((pm3d.clip == PM3D_CLIP_Z)
1584 && (pm3d_plot_at != PM3D_AT_BASE && pm3d_plot_at != PM3D_AT_TOP)) {
1585 int new = clip_filled_polygon( corners, clipcorners, nv );
1586 if (new < 0) /* All vertices out of range */
1587 return;
1588 if (new > 0) { /* Some got clipped */
1589 nv = new;
1590 corners = clipcorners;
1591 }
1592 }
1593
1594 for (i = 0; i < nv; i++) {
1595 map3d_xy_double(corners[i].x, corners[i].y, corners[i].z, &x, &y);
1596 icorners[i].x = x;
1597 icorners[i].y = y;
1598 }
1599
1600 /* Clip to x/y only in 2D projection. */
1601 if (splot_map && (pm3d.clip == PM3D_CLIP_Z)) {
1602 for (i=0; i<nv; i++)
1603 ocorners[i] = icorners[i];
1604 clip_polygon( ocorners, icorners, nv, &nv );
1605 }
1606
1607 if (fillstyle)
1608 icorners[0].style = fillstyle;
1609 else if (default_fillstyle.fillstyle == FS_EMPTY)
1610 icorners[0].style = FS_OPAQUE;
1611 else
1612 icorners[0].style = style_from_fill(&default_fillstyle);
1613
1614 term->filled_polygon(nv, icorners);
1615
1616 if (pm3d.border.l_type != LT_NODRAW) {
1617 /* LT_DEFAULT means draw border in current color */
1618 /* FIXME: currently there is no obvious way to set LT_DEFAULT */
1619 if (pm3d.border.l_type != LT_DEFAULT)
1620 term_apply_lp_properties(&pm3d.border);
1621
1622 term->move(icorners[0].x, icorners[0].y);
1623 for (i = nv-1; i >= 0; i--) {
1624 term->vector(icorners[i].x, icorners[i].y);
1625 }
1626 }
1627 }
1628
1629 /*
1630 * Clip existing filled polygon again zmax or zmin.
1631 * We already know this is a convex polygon with at least one vertex in range.
1632 * The clipped polygon may have as few as 3 vertices or as many as n+2.
1633 * Returns the new number of vertices after clipping.
1634 */
1635 int
clip_filled_polygon(gpdPoint * inpts,gpdPoint * outpts,int nv)1636 clip_filled_polygon( gpdPoint *inpts, gpdPoint *outpts, int nv )
1637 {
1638 int current = 0; /* The vertex we are now considering */
1639 int next = 0; /* The next vertex */
1640 int nvo = 0; /* Number of vertices after clipping */
1641 int noutrange = 0; /* Number of out-range points */
1642 int nover = 0;
1643 int nunder = 0;
1644 double fraction;
1645 double zmin = axis_array[FIRST_Z_AXIS].min;
1646 double zmax = axis_array[FIRST_Z_AXIS].max;
1647
1648 /* classify inrange/outrange vertices */
1649 static int *outrange = NULL;
1650 static int maxvert = 0;
1651 if (nv > maxvert) {
1652 maxvert = nv;
1653 outrange = gp_realloc(outrange, maxvert * sizeof(int), NULL);
1654 }
1655 for (current = 0; current < nv; current++) {
1656 if (inrange( inpts[current].z, zmin, zmax )) {
1657 outrange[current] = 0;
1658 } else if (inpts[current].z > zmax) {
1659 outrange[current] = 1;
1660 noutrange++;
1661 nover++;
1662 } else if (inpts[current].z < zmin) {
1663 outrange[current] = -1;
1664 noutrange++;
1665 nunder++;
1666 }
1667 }
1668
1669 /* If all are in range, nothing to do */
1670 if (noutrange == 0)
1671 return 0;
1672
1673 /* If all vertices are out of range on the same side, draw nothing */
1674 if (nunder == nv || nover == nv)
1675 return -1;
1676
1677 for (current = 0; current < nv; current++) {
1678 next = (current+1) % nv;
1679
1680 /* Current point is out of range */
1681 if (outrange[current]) {
1682 if (!outrange[next]) {
1683 /* Current point is out-range but next point is in-range.
1684 * Clip line segment from current-to-next and store as new vertex.
1685 */
1686 fraction = ((inpts[current].z >= zmax ? zmax : zmin) - inpts[next].z)
1687 / (inpts[current].z - inpts[next].z);
1688 outpts[nvo].x = inpts[next].x + fraction * (inpts[current].x - inpts[next].x);
1689 outpts[nvo].y = inpts[next].y + fraction * (inpts[current].y - inpts[next].y);
1690 outpts[nvo].z = inpts[current].z >= zmax ? zmax : zmin;
1691 nvo++;
1692 }
1693 else if ( /* clip_lines2 && */ (outrange[current] * outrange[next] < 0)) {
1694 /* Current point and next point are out of range on opposite
1695 * sides of the z range. Clip both ends of the segment.
1696 */
1697 fraction = ((inpts[current].z >= zmax ? zmax : zmin) - inpts[next].z)
1698 / (inpts[current].z - inpts[next].z);
1699 outpts[nvo].x = inpts[next].x + fraction * (inpts[current].x - inpts[next].x);
1700 outpts[nvo].y = inpts[next].y + fraction * (inpts[current].y - inpts[next].y);
1701 outpts[nvo].z = inpts[current].z >= zmax ? zmax : zmin;
1702 nvo++;
1703 fraction = ((inpts[next].z >= zmax ? zmax : zmin) - outpts[nvo-1].z)
1704 / (inpts[next].z - outpts[nvo-1].z);
1705 outpts[nvo].x = outpts[nvo-1].x + fraction * (inpts[next].x - outpts[nvo-1].x);
1706 outpts[nvo].y = outpts[nvo-1].y + fraction * (inpts[next].y - outpts[nvo-1].y);
1707 outpts[nvo].z = inpts[next].z >= zmax ? zmax : zmin;
1708 nvo++;
1709 }
1710
1711 /* Current point is in range */
1712 } else {
1713 outpts[nvo++] = inpts[current];
1714 if (outrange[next]) {
1715 /* Current point is in-range but next point is out-range.
1716 * Clip line segment from current-to-next and store as new vertex.
1717 */
1718 fraction = ((inpts[next].z >= zmax ? zmax : zmin) - inpts[current].z)
1719 / (inpts[next].z - inpts[current].z);
1720 outpts[nvo].x = inpts[current].x + fraction * (inpts[next].x - inpts[current].x);
1721 outpts[nvo].y = inpts[current].y + fraction * (inpts[next].y - inpts[current].y);
1722 outpts[nvo].z = inpts[next].z >= zmax ? zmax : zmin;
1723 nvo++;
1724 }
1725 }
1726 }
1727
1728 /* this is not supposed to happen, but ... */
1729 if (nvo <= 1)
1730 return -1;
1731
1732 return nvo;
1733 }
1734
1735
1736 /* EXPERIMENTAL
1737 * returns 1 for top of pm3d surface towards the viewer
1738 * -1 for bottom of pm3d surface towards the viewer
1739 * NB: the ordering of the quadrangle vertices depends on the scan direction.
1740 * In the case of depth ordering, the user does not have good control
1741 * over this.
1742 */
1743 int
pm3d_side(struct coordinate * p0,struct coordinate * p1,struct coordinate * p2)1744 pm3d_side( struct coordinate *p0, struct coordinate *p1, struct coordinate *p2)
1745 {
1746 struct vertex v[3];
1747 double u0, u1, v0, v1;
1748
1749 /* Apply current view rotation to corners of this quadrangle */
1750 map3d_xyz(p0->x, p0->y, p0->z, &v[0]);
1751 map3d_xyz(p1->x, p1->y, p1->z, &v[1]);
1752 map3d_xyz(p2->x, p2->y, p2->z, &v[2]);
1753
1754 /* projection of two adjacent edges */
1755 u0 = v[1].x - v[0].x;
1756 u1 = v[1].y - v[0].y;
1757 v0 = v[2].x - v[0].x;
1758 v1 = v[2].y - v[0].y;
1759
1760 /* cross-product */
1761 return sgn(u0*v1 - u1*v0);
1762 }
1763
1764 /*
1765 * Returns a pointer into the list of polygon vertices.
1766 * If necessary reallocates the entire list to ensure there is enough
1767 * room for the requested number of vertices.
1768 */
get_polygon(int size)1769 static gpdPoint *get_polygon( int size )
1770 {
1771 /* Check size of current list, allocate more space if needed */
1772 if (next_polygon + size >= polygonlistsize) {
1773 polygonlistsize = 2*polygonlistsize + size;
1774 polygonlist = gp_realloc(polygonlist, polygonlistsize * sizeof(gpdPoint), NULL);
1775 }
1776
1777 current_polygon = next_polygon;
1778 next_polygon = current_polygon + size;
1779 return &polygonlist[current_polygon];
1780 }
1781
1782 void
free_polygonlist()1783 free_polygonlist()
1784 {
1785 free(polygonlist);
1786 polygonlist = NULL;
1787 current_polygon = 0;
1788 next_polygon = 0;
1789 polygonlistsize = 0;
1790 }
1791
1792 void
pm3d_reset_after_error()1793 pm3d_reset_after_error()
1794 {
1795 pm3d_plot_at = 0;
1796 free_polygonlist();
1797 }
1798