1 /* $Id: plshade.c,v 1.3 2007/05/08 09:09:37 rice Exp $
2
3 Functions to shade regions on the basis of value.
4 Can be used to shade contour plots or alone.
5 Copyright 1993 Wesley Ebisuzaki
6
7 Copyright (C) 2004 Alan W. Irwin
8
9 This file is part of PLplot.
10
11 PLplot is free software; you can redistribute it and/or modify
12 it under the terms of the GNU General Library Public License as published
13 by the Free Software Foundation; either version 2 of the License, or
14 (at your option) any later version.
15
16 PLplot is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 GNU Library General Public License for more details.
20
21 You should have received a copy of the GNU Library General Public License
22 along with PLplot; if not, write to the Free Software
23 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24
25 */
26
27 /*----------------------------------------------------------------------*\
28 * Call syntax for plshade():
29 *
30 * void plshade(PLFLT *a, PLINT nx, PLINT ny, char *defined,
31 * PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
32 * PLFLT shade_min, PLFLT shade_max,
33 * PLINT sh_color, PLINT sh_width, PLINT min_color, PLINT min_width,
34 * PLINT max_color, PLINT max_width, void (*fill)(), PLINT
35 * rectangular, ...)
36 *
37 * arguments:
38 *
39 * PLFLT &(a[0][0])
40 *
41 * Contains array to be plotted. The array must have been declared as
42 * PLFLT a[nx][ny]. See following note on fortran-style arrays.
43 *
44 * PLINT nx, ny
45 *
46 * Dimension of array "a".
47 *
48 * char &(defined[0][0])
49 *
50 * Contains array of flags, 1 = data is valid, 0 = data is not valid.
51 * This array determines which sections of the data is to be plotted.
52 * This argument can be NULL if all the values are valid. Must have been
53 * declared as char defined[nx][ny].
54 *
55 * PLFLT xmin, xmax, ymin, ymax
56 *
57 * Defines the "grid" coordinates. The data a[0][0] has a position of
58 * (xmin,ymin).
59 *
60 * void (*mapform)()
61 *
62 * Transformation from `grid' coordinates to world coordinates. This
63 * pointer to a function can be NULL in which case the grid coordinates
64 * are the same as the world coordinates.
65 *
66 * PLFLT shade_min, shade_max
67 *
68 * Defines the interval to be shaded. If shade_max <= shade_min, plshade
69 * does nothing.
70 *
71 * PLINT sh_cmap, PLFLT sh_color, PLINT sh_width
72 *
73 * Defines color map, color map index, and width used by the fill pattern.
74 *
75 * PLINT min_color, min_width, max_color, max_width
76 *
77 * Defines pen color, width used by the boundary of shaded region. The min
78 * values are used for the shade_min boundary, and the max values are used
79 * on the shade_max boundary. Set color and width to zero for no plotted
80 * boundaries.
81 *
82 * void (*fill)()
83 *
84 * Routine used to fill the region. Use plfill. Future version of plplot
85 * may have other fill routines.
86 *
87 * PLINT rectangular
88 *
89 * Flag. Set to 1 if rectangles map to rectangles after (*mapform)() else
90 * set to zero. If rectangular is set to 1, plshade tries to save time by
91 * filling large rectangles. This optimization fails if (*mapform)()
92 * distorts the shape of rectangles. For example a plot in polor
93 * coordinates has to have rectangular set to zero.
94 *
95 * Example mapform's:
96 *
97 * Grid to world coordinate transformation.
98 * This example goes from a r-theta to x-y for a polar plot.
99 *
100 * void mapform(PLINT n, PLFLT *x, PLFLT *y) {
101 * int i;
102 * double r, theta;
103 * for (i = 0; i < n; i++) {
104 * r = x[i];
105 * theta = y[i];
106 * x[i] = r*cos(theta);
107 * y[i] = r*sin(theta);
108 * }
109 * }
110 *
111 * Grid was in cm, convert to world coordinates in inches.
112 * Expands in x direction.
113 *
114 * void mapform(PLINT n, PLFLT *x, PLFLT *y) {
115 * int i;
116 * for (i = 0; i < n; i++) {
117 * x[i] = (1.0 / 2.5) * x[i];
118 * y[i] = (1.0 / 2.5) * y[i];
119 * }
120 * }
121 *
122 \*----------------------------------------------------------------------*/
123
124 #include "plplotP.h"
125 #include <float.h>
126
127 #define MISSING_MIN_DEF (PLFLT) 1.0
128 #define MISSING_MAX_DEF (PLFLT) -1.0
129
130
131 #define NEG 1
132 #define POS 8
133 #define OK 0
134 #define UNDEF 64
135
136 #define linear(val1, val2, level) ((level - val1) / (val2 - val1))
137
138 /* Global variables */
139
140 static PLFLT sh_max, sh_min;
141 static int min_points, max_points, n_point;
142 static int min_pts[4], max_pts[4];
143 static PLINT pen_col_min, pen_col_max;
144 static PLINT pen_wd_min, pen_wd_max;
145 static PLFLT int_val;
146
147 /* Function prototypes */
148
149 static void
150 set_cond(register int *cond, register PLFLT *a, register PLINT n);
151
152 static int
153 find_interval(PLFLT a0, PLFLT a1, PLINT c0, PLINT c1, PLFLT *x);
154
155 static void
156 selected_polygon(void (*fill) (PLINT, PLFLT *, PLFLT *),
157 PLINT (*defined) (PLFLT, PLFLT),
158 PLFLT *x, PLFLT *y, PLINT v1, PLINT v2, PLINT v3, PLINT v4);
159
160 static void
161 exfill(void (*fill) (PLINT, PLFLT *, PLFLT *),
162 PLINT (*defined) (PLFLT, PLFLT),
163 int n, PLFLT *x, PLFLT *y);
164
165 static void
166 big_recl(int *cond_code, register int ny, int dx, int dy,
167 int *ix, int *iy);
168
169 static void
170 draw_boundary(PLINT slope, PLFLT *x, PLFLT *y);
171
172 static PLINT
173 plctest(PLFLT *x, PLFLT level);
174
175 static PLINT
176 plctestez(PLFLT *a, PLINT nx, PLINT ny, PLINT ix,
177 PLINT iy, PLFLT level);
178
179 static void
180 plshade_int(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
181 PLPointer f2eval_data,
182 PLFLT (*c2eval) (PLINT, PLINT, PLPointer),
183 PLPointer c2eval_data,
184 PLINT (*defined) (PLFLT, PLFLT),
185 PLFLT missing_min, PLFLT missing_max,
186 PLINT nx, PLINT ny,
187 PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
188 PLFLT shade_min, PLFLT shade_max,
189 PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
190 PLINT min_color, PLINT min_width,
191 PLINT max_color, PLINT max_width,
192 void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
193 void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
194 PLPointer pltr_data);
195
196 /*----------------------------------------------------------------------*\
197 * plshades()
198 *
199 * Shade regions via a series of calls to plshade.
200 * All arguments are the same as plshade except the following:
201 * clevel is a pointer to an array of values representing
202 * the shade edge values, nlevel-1 is
203 * the number of different shades, (nlevel is the number of shade edges),
204 * fill_width is the pattern fill width, and cont_color and cont_width
205 * are the color and width of the contour drawn at each shade edge.
206 * (if cont_color <= 0 or cont_width <=0, no such contours are drawn).
207 \*----------------------------------------------------------------------*/
208
c_plshades(PLFLT ** a,PLINT nx,PLINT ny,PLINT (* defined)(PLFLT,PLFLT),PLFLT xmin,PLFLT xmax,PLFLT ymin,PLFLT ymax,PLFLT * clevel,PLINT nlevel,PLINT fill_width,PLINT cont_color,PLINT cont_width,void (* fill)(PLINT,PLFLT *,PLFLT *),PLINT rectangular,void (* pltr)(PLFLT,PLFLT,PLFLT *,PLFLT *,PLPointer),PLPointer pltr_data)209 void c_plshades( PLFLT **a, PLINT nx, PLINT ny, PLINT (*defined) (PLFLT, PLFLT),
210 PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
211 PLFLT *clevel, PLINT nlevel, PLINT fill_width,
212 PLINT cont_color, PLINT cont_width,
213 void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
214 void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
215 PLPointer pltr_data )
216 {
217 PLFLT shade_min, shade_max, shade_color;
218 PLINT i, init_color, init_width;
219
220 for (i = 0; i < nlevel-1; i++) {
221 shade_min = clevel[i];
222 shade_max = clevel[i+1];
223 shade_color = i / (PLFLT) (nlevel-2);
224 /* The constants in order mean
225 * (1) color map1,
226 * (0, 0, 0, 0) all edge effects will be done with plcont rather
227 * than the normal plshade drawing which gets partially blocked
228 * when sequential shading is done as in the present case */
229
230 plshade(a, nx, ny, defined, xmin, xmax, ymin, ymax,
231 shade_min, shade_max,
232 1, shade_color, fill_width,
233 0, 0, 0, 0,
234 fill, rectangular, pltr, pltr_data);
235 }
236 if(cont_color > 0 && cont_width > 0) {
237 init_color = plsc->icol0;
238 init_width = plsc->width;
239 plcol0(cont_color);
240 plwid(cont_width);
241 if(pltr && pltr_data) {
242 plcont(a, nx, ny, 1, nx, 1, ny, clevel, nlevel, pltr, pltr_data);
243 }
244 else {
245 /* For this case use the same interpretation that occurs internally
246 * for plshade. That is set up x and y grids that map from the
247 * index ranges to xmin, xmax, ymin, ymax, and use those grids
248 * for the plcont call.
249 */
250 PLcGrid cgrid1;
251 PLFLT *x, *y;
252 cgrid1.nx = nx;
253 cgrid1.ny = ny;
254 x = (PLFLT *) malloc(nx * sizeof(PLFLT));
255 if (x == NULL)
256 plexit("plshades: Out of memory for x");
257 cgrid1.xg = x;
258 for (i = 0; i < nx; i++)
259 cgrid1.xg[i] = xmin + (xmax - xmin)*(float)i/(float)(nx-1);
260 y = (PLFLT *) malloc(ny * sizeof(PLFLT));
261 if (y == NULL)
262 plexit("plshades: Out of memory for y");
263 cgrid1.yg = y;
264 for (i = 0; i < ny; i++)
265 cgrid1.yg[i] = ymin + (ymax - ymin)*(float)i/(float)(ny-1);
266 plcont(a, nx, ny, 1, nx, 1, ny, clevel, nlevel,
267 pltr1, (void *) &cgrid1);
268 free(x);
269 free(y);
270 }
271 plcol0(init_color);
272 plwid(init_width);
273 }
274 }
275
276 /*----------------------------------------------------------------------*\
277 * plshade()
278 *
279 * Shade region.
280 * This interface to plfshade() assumes the 2d function array is passed
281 * via a (PLFLT **), and is column-dominant (normal C ordering).
282 \*----------------------------------------------------------------------*/
283
c_plshade(PLFLT ** a,PLINT nx,PLINT ny,PLINT (* defined)(PLFLT,PLFLT),PLFLT xmin,PLFLT xmax,PLFLT ymin,PLFLT ymax,PLFLT shade_min,PLFLT shade_max,PLINT sh_cmap,PLFLT sh_color,PLINT sh_width,PLINT min_color,PLINT min_width,PLINT max_color,PLINT max_width,void (* fill)(PLINT,PLFLT *,PLFLT *),PLINT rectangular,void (* pltr)(PLFLT,PLFLT,PLFLT *,PLFLT *,PLPointer),PLPointer pltr_data)284 void c_plshade( PLFLT **a, PLINT nx, PLINT ny, PLINT (*defined) (PLFLT, PLFLT),
285 PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
286 PLFLT shade_min, PLFLT shade_max,
287 PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
288 PLINT min_color, PLINT min_width,
289 PLINT max_color, PLINT max_width,
290 void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
291 void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
292 PLPointer pltr_data )
293 {
294 PLfGrid2 grid;
295
296 grid.f = a;
297 grid.nx = nx;
298 grid.ny = ny;
299
300 plshade_int( plf2eval2, (PLPointer) &grid,
301 NULL, NULL,
302 /* plc2eval, (PLPointer) &cgrid,*/
303 defined, MISSING_MIN_DEF, MISSING_MAX_DEF, nx, ny, xmin,
304 xmax, ymin, ymax, shade_min, shade_max,
305 sh_cmap, sh_color, sh_width,
306 min_color, min_width, max_color, max_width,
307 fill, rectangular, pltr, pltr_data );
308 }
309
310 /*----------------------------------------------------------------------*\
311 * plshade1()
312 *
313 * Shade region.
314 * This interface to plfshade() assumes the 2d function array is passed
315 * via a (PLFLT *), and is column-dominant (normal C ordering).
316 \*----------------------------------------------------------------------*/
317
c_plshade1(PLFLT * a,PLINT nx,PLINT ny,PLINT (* defined)(PLFLT,PLFLT),PLFLT xmin,PLFLT xmax,PLFLT ymin,PLFLT ymax,PLFLT shade_min,PLFLT shade_max,PLINT sh_cmap,PLFLT sh_color,PLINT sh_width,PLINT min_color,PLINT min_width,PLINT max_color,PLINT max_width,void (* fill)(PLINT,PLFLT *,PLFLT *),PLINT rectangular,void (* pltr)(PLFLT,PLFLT,PLFLT *,PLFLT *,PLPointer),PLPointer pltr_data)318 void c_plshade1( PLFLT *a, PLINT nx, PLINT ny, PLINT (*defined) (PLFLT, PLFLT),
319 PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
320 PLFLT shade_min, PLFLT shade_max,
321 PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
322 PLINT min_color, PLINT min_width,
323 PLINT max_color, PLINT max_width,
324 void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
325 void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
326 PLPointer pltr_data )
327 {
328 PLfGrid grid;
329
330 grid.f = a;
331 grid.nx = nx;
332 grid.ny = ny;
333
334 plshade_int( plf2eval, (PLPointer) &grid,
335 NULL, NULL,
336 /* plc2eval, (PLPointer) &cgrid,*/
337 defined, MISSING_MIN_DEF, MISSING_MAX_DEF, nx, ny, xmin,
338 xmax, ymin, ymax, shade_min, shade_max,
339 sh_cmap, sh_color, sh_width,
340 min_color, min_width, max_color, max_width,
341 fill, rectangular, pltr, pltr_data );
342 }
343
344 /*----------------------------------------------------------------------*\
345 * plfshade()
346 *
347 * Shade region.
348 * Array values are determined by the input function and the passed data.
349 \*----------------------------------------------------------------------*/
350
351 void
plfshade(PLFLT (* f2eval)(PLINT,PLINT,PLPointer),PLPointer f2eval_data,PLFLT (* c2eval)(PLINT,PLINT,PLPointer),PLPointer c2eval_data,PLINT nx,PLINT ny,PLFLT xmin,PLFLT xmax,PLFLT ymin,PLFLT ymax,PLFLT shade_min,PLFLT shade_max,PLINT sh_cmap,PLFLT sh_color,PLINT sh_width,PLINT min_color,PLINT min_width,PLINT max_color,PLINT max_width,void (* fill)(PLINT,PLFLT *,PLFLT *),PLINT rectangular,void (* pltr)(PLFLT,PLFLT,PLFLT *,PLFLT *,PLPointer),PLPointer pltr_data)352 plfshade(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
353 PLPointer f2eval_data,
354 PLFLT (*c2eval) (PLINT, PLINT, PLPointer),
355 PLPointer c2eval_data,
356 PLINT nx, PLINT ny,
357 PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
358 PLFLT shade_min, PLFLT shade_max,
359 PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
360 PLINT min_color, PLINT min_width,
361 PLINT max_color, PLINT max_width,
362 void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
363 void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
364 PLPointer pltr_data)
365 {
366 plshade_int(f2eval, f2eval_data, c2eval, c2eval_data,
367 NULL, MISSING_MIN_DEF, MISSING_MAX_DEF,
368 nx, ny, xmin, xmax, ymin, ymax,
369 shade_min, shade_max, sh_cmap, sh_color, sh_width,
370 min_color, min_width, max_color, max_width,
371 fill, rectangular, pltr, pltr_data);
372 }
373
374
375 /*----------------------------------------------------------------------*\
376 * plshade_int()
377 *
378 * Shade region -- this routine does all the work
379 *
380 * This routine is internal so the arguments can and will change.
381 * To retain some compatibility between versions, you must go through
382 * some stub routine!
383 *
384 * 4/95
385 *
386 * new: missing_min, missing_max
387 *
388 * if data <= missing_max and data >= missing_min
389 * then the data will beconsidered to be missing
390 * this allows 2nd way to set undefined points (good for ftn)
391 * if missing feature is not used, set missing_max < missing_min
392 *
393 * parameters:
394 *
395 * f2eval, f2eval_data: data to plot
396 * c2eval, c2eval_data: defined mask (not implimented)
397 * defined: defined mask (old API - implimented)
398 * missing_min, missing_max: yet another way to set data to undefined
399 * nx, ny: array dimensions
400 * xmin, xmax, ymin, ymax: grid coordinates
401 * shade_min, shade_max: shade region with values between ...
402 * sh_cmap, sh_color, sh_width: shading parameters, width is only for hatching
403 * min_color, min_width: line parameters for boundary (minimum)
404 * max_color, max_width: line parameters for boundary (maximum)
405 * set min_width == 0 and max_width == 0 for no contours
406 * fill: fill function, set to NULL for no shading (contour plot)
407 * rectangular: flag set to 1 if pltr() maps rectangles to rectangles
408 * this helps optimize the plotting
409 * pltr: function to map from grid to plot coordinates
410 *
411 *
412 \*----------------------------------------------------------------------*/
413
414 static void
plshade_int(PLFLT (* f2eval)(PLINT,PLINT,PLPointer),PLPointer f2eval_data,PLFLT (* c2eval)(PLINT,PLINT,PLPointer),PLPointer c2eval_data,PLINT (* defined)(PLFLT,PLFLT),PLFLT missing_min,PLFLT missing_max,PLINT nx,PLINT ny,PLFLT xmin,PLFLT xmax,PLFLT ymin,PLFLT ymax,PLFLT shade_min,PLFLT shade_max,PLINT sh_cmap,PLFLT sh_color,PLINT sh_width,PLINT min_color,PLINT min_width,PLINT max_color,PLINT max_width,void (* fill)(PLINT,PLFLT *,PLFLT *),PLINT rectangular,void (* pltr)(PLFLT,PLFLT,PLFLT *,PLFLT *,PLPointer),PLPointer pltr_data)415 plshade_int(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
416 PLPointer f2eval_data,
417 PLFLT (*c2eval) (PLINT, PLINT, PLPointer),
418 PLPointer c2eval_data,
419 PLINT (*defined) (PLFLT, PLFLT),
420 PLFLT missing_min, PLFLT missing_max,
421 PLINT nx, PLINT ny,
422 PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
423 PLFLT shade_min, PLFLT shade_max,
424 PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
425 PLINT min_color, PLINT min_width,
426 PLINT max_color, PLINT max_width,
427 void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
428 void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
429 PLPointer pltr_data)
430 {
431
432 PLINT init_width, n, slope = 0, ix, iy;
433 int count, i, j, nxny;
434 PLFLT *a, *a0, *a1, dx, dy;
435 PLFLT x[8], y[8], xp[2], tx, ty;
436 int *c, *c0, *c1;
437
438 (void) c2eval; /* pmr: make it used */
439 (void) c2eval_data;
440 (void) missing_min;
441 (void) missing_max;
442
443 if (plsc->level < 3) {
444 plabort("plfshade: window must be set up first");
445 return;
446 }
447
448 if (nx <= 0 || ny <= 0) {
449 plabort("plfshade: nx and ny must be positive");
450 return;
451 }
452
453 if (shade_min >= shade_max) {
454 plabort("plfshade: shade_max must exceed shade_min");
455 return;
456 }
457
458 if (pltr == NULL || pltr_data == NULL)
459 rectangular = 1;
460
461 int_val = shade_max - shade_min;
462 init_width = plsc->width;
463
464 pen_col_min = min_color;
465 pen_col_max = max_color;
466
467 pen_wd_min = min_width;
468 pen_wd_max = max_width;
469
470 plstyl((PLINT) 0, NULL, NULL);
471 plwid(sh_width);
472 if (fill != NULL) {
473 switch (sh_cmap) {
474 case 0:
475 plcol0((PLINT) sh_color);
476 break;
477 case 1:
478 plcol1(sh_color);
479 break;
480 default:
481 plabort("plfshade: invalid color map selection");
482 return;
483 }
484 }
485 /* alloc space for value array, and initialize */
486 /* This is only a temporary kludge */
487 nxny = nx * ny;
488 if ((a = (PLFLT *) malloc(nxny * sizeof(PLFLT))) == NULL) {
489 plabort("plfshade: unable to allocate memory for value array");
490 return;
491 }
492
493 for (ix = 0; ix < nx; ix++)
494 for (iy = 0; iy < ny; iy++)
495 a[iy + ix*ny] = f2eval(ix, iy, f2eval_data);
496
497 /* alloc space for condition codes */
498
499 if ((c = (int *) malloc(nxny * sizeof(int))) == NULL) {
500 plabort("plfshade: unable to allocate memory for condition codes");
501 free(a);
502 return;
503 }
504
505 sh_min = shade_min;
506 sh_max = shade_max;
507
508 set_cond(c, a, nxny);
509 dx = (xmax - xmin) / (nx - 1);
510 dy = (ymax - ymin) / (ny - 1);
511 a0 = a;
512 a1 = a + ny;
513 c0 = c;
514 c1 = c + ny;
515
516 for (ix = 0; ix < nx - 1; ix++) {
517
518 for (iy = 0; iy < ny - 1; iy++) {
519
520 count = c0[iy] + c0[iy + 1] + c1[iy] + c1[iy + 1];
521
522 /* No filling needs to be done for these cases */
523
524 if (count >= UNDEF)
525 continue;
526 if (count == 4 * POS)
527 continue;
528 if (count == 4 * NEG)
529 continue;
530
531 /* Entire rectangle can be filled */
532
533 if (count == 4 * OK) {
534 /* find biggest rectangle that fits */
535 if (rectangular) {
536 big_recl(c0 + iy, ny, nx - ix, ny - iy, &i, &j);
537 }
538 else {
539 i = j = 1;
540 }
541 x[0] = x[1] = ix;
542 x[2] = x[3] = ix+i;
543 y[0] = y[3] = iy;
544 y[1] = y[2] = iy+j;
545
546 if (pltr && pltr_data) {
547 for (i = 0; i < 4; i++) {
548 (*pltr) (x[i], y[i], &tx, &ty, pltr_data);
549 x[i] = tx;
550 y[i] = ty;
551 }
552 }
553 else {
554 for (i = 0; i < 4; i++) {
555 x[i] = xmin + x[i]*dx;
556 y[i] = ymin + y[i]*dy;
557 }
558 }
559 if (fill != NULL)
560 exfill (fill, defined, (PLINT) 4, x, y);
561 iy += j - 1;
562 continue;
563 }
564
565 /* Only part of rectangle can be filled */
566
567 n_point = min_points = max_points = 0;
568 n = find_interval(a0[iy], a0[iy + 1], c0[iy], c0[iy + 1], xp);
569 for (j = 0; j < n; j++) {
570 x[j] = ix;
571 y[j] = iy + xp[j];
572 }
573
574 i = find_interval(a0[iy + 1], a1[iy + 1],
575 c0[iy + 1], c1[iy + 1], xp);
576
577 for (j = 0; j < i; j++) {
578 x[j + n] = ix + xp[j];
579 y[j + n] = iy + 1;
580 }
581 n += i;
582
583 i = find_interval(a1[iy + 1], a1[iy], c1[iy + 1], c1[iy], xp);
584 for (j = 0; j < i; j++) {
585 x[n + j] = ix + 1;
586 y[n + j] = iy + 1 - xp[j];
587 }
588 n += i;
589
590 i = find_interval(a1[iy], a0[iy], c1[iy], c0[iy], xp);
591 for (j = 0; j < i; j++) {
592 x[n + j] = ix + 1 - xp[j];
593 y[n + j] = iy;
594 }
595 n += i;
596
597 if (pltr && pltr_data) {
598 for (i = 0; i < n; i++) {
599 (*pltr) (x[i], y[i], &tx, &ty, pltr_data);
600 x[i] = tx;
601 y[i] = ty;
602 }
603 }
604 else {
605 for (i = 0; i < n; i++) {
606 x[i] = xmin + x[i]*dx;
607 y[i] = ymin + y[i]*dy;
608 }
609 }
610
611 if (min_points == 4)
612 slope = plctestez(a, nx, ny, ix, iy, shade_min);
613 if (max_points == 4)
614 slope = plctestez(a, nx, ny, ix, iy, shade_max);
615
616 /* n = number of end of line segments */
617 /* min_points = number times shade_min meets edge */
618 /* max_points = number times shade_max meets edge */
619
620 /* special cases: check number of times a contour is in a box */
621
622 switch ((min_points << 3) + max_points) {
623 case 000:
624 case 020:
625 case 002:
626 case 022:
627 if (fill != NULL && n > 0)
628 exfill (fill, defined, n, x, y);
629 break;
630 case 040: /* 2 contour lines in box */
631 case 004:
632 if (n != 6)
633 fprintf(stderr, "plfshade err n=%d !6", (int) n);
634 if (slope == 1 && c0[iy] == OK) {
635 if (fill != NULL)
636 exfill (fill, defined, n, x, y);
637 }
638 else if (slope == 1) {
639 selected_polygon(fill, defined, x, y, 0, 1, 2, -1);
640 selected_polygon(fill, defined, x, y, 3, 4, 5, -1);
641 }
642 else if (c0[iy + 1] == OK) {
643 if (fill != NULL)
644 exfill (fill, defined, n, x, y);
645 }
646 else {
647 selected_polygon(fill, defined, x, y, 0, 1, 5, -1);
648 selected_polygon(fill, defined, x, y, 2, 3, 4, -1);
649 }
650 break;
651 case 044:
652 if (n != 8)
653 fprintf(stderr, "plfshade err n=%d !8", (int) n);
654 if (slope == 1) {
655 selected_polygon(fill, defined, x, y, 0, 1, 2, 3);
656 selected_polygon(fill, defined, x, y, 4, 5, 6, 7);
657 }
658 else {
659 selected_polygon(fill, defined, x, y, 0, 1, 6, 7);
660 selected_polygon(fill, defined, x, y, 2, 3, 4, 5);
661 }
662 break;
663 case 024:
664 case 042:
665 /* 3 contours */
666 if (n != 7)
667 fprintf(stderr, "plfshade err n=%d !7", (int) n);
668
669 if ((c0[iy] == OK || c1[iy+1] == OK) && slope == 1) {
670 if (fill != NULL)
671 exfill (fill, defined, n, x, y);
672 }
673 else if ((c0[iy+1] == OK || c1[iy] == OK) && slope == 0) {
674 if (fill !=NULL)
675 exfill (fill, defined, n, x, y);
676 }
677
678 else if (c0[iy] == OK) {
679 selected_polygon(fill, defined, x, y, 0, 1, 6, -1);
680 selected_polygon(fill, defined, x, y, 2, 3, 4, 5);
681 }
682 else if (c0[iy+1] == OK) {
683 selected_polygon(fill, defined, x, y, 0, 1, 2, -1);
684 selected_polygon(fill, defined, x, y, 3, 4, 5, 6);
685 }
686 else if (c1[iy+1] == OK) {
687 selected_polygon(fill, defined, x, y, 0, 1, 5, 6);
688 selected_polygon(fill, defined, x, y, 2, 3, 4, -1);
689 }
690 else if (c1[iy] == OK) {
691 selected_polygon(fill, defined, x, y, 0, 1, 2, 3);
692 selected_polygon(fill, defined, x, y, 4, 5, 6, -1);
693 }
694 else {
695 fprintf(stderr, "plfshade err logic case 024:042\n");
696 }
697 break;
698 default:
699 fprintf(stderr, "prog err switch\n");
700 break;
701 }
702 draw_boundary(slope, x, y);
703
704 if (fill != NULL) {
705 plwid(sh_width);
706 if (sh_cmap == 0) plcol0((PLINT) sh_color);
707 else if (sh_cmap == 1) plcol1(sh_color);
708 }
709 }
710
711 a0 = a1;
712 c0 = c1;
713 a1 += ny;
714 c1 += ny;
715 }
716
717 free(c);
718 free(a);
719 plwid(init_width);
720 }
721
722 /*----------------------------------------------------------------------*\
723 * set_cond()
724 *
725 * Fills out condition code array.
726 \*----------------------------------------------------------------------*/
727
728 static void
set_cond(register int * cond,register PLFLT * a,register PLINT n)729 set_cond(register int *cond, register PLFLT *a, register PLINT n)
730 {
731 while (n--) {
732 if (*a < sh_min)
733 *cond++ = NEG;
734 else if (*a > sh_max)
735 *cond++ = POS;
736 else
737 *cond++ = OK;
738 a++;
739 }
740 }
741
742 /*----------------------------------------------------------------------*\
743 * find_interval()
744 *
745 * Two points x(0) = a0, (condition code c0) x(1) = a1, (condition code c1)
746 * return interval on the line that are shaded
747 *
748 * returns 0 : no points to be shaded 1 : x[0] <= x < 1 is the interval 2 :
749 * x[0] <= x <= x[1] < 1 interval to be shaded n_point, max_points,
750 * min_points are incremented location of min/max_points are stored
751 \*----------------------------------------------------------------------*/
752
753 static int
find_interval(PLFLT a0,PLFLT a1,PLINT c0,PLINT c1,PLFLT * x)754 find_interval(PLFLT a0, PLFLT a1, PLINT c0, PLINT c1, PLFLT *x)
755 {
756 register int n;
757
758 n = 0;
759 if (c0 == OK) {
760 x[n++] = 0.0;
761 n_point++;
762 }
763 if (c0 == c1)
764 return n;
765
766 if (c0 == NEG || c1 == POS) {
767 if (c0 == NEG) {
768 x[n++] = linear(a0, a1, sh_min);
769 min_pts[min_points++] = n_point++;
770 }
771 if (c1 == POS) {
772 x[n++] = linear(a0, a1, sh_max);
773 max_pts[max_points++] = n_point++;
774 }
775 }
776 if (c0 == POS || c1 == NEG) {
777 if (c0 == POS) {
778 x[n++] = linear(a0, a1, sh_max);
779 max_pts[max_points++] = n_point++;
780 }
781 if (c1 == NEG) {
782 x[n++] = linear(a0, a1, sh_min);
783 min_pts[min_points++] = n_point++;
784 }
785 }
786 return n;
787 }
788
789 /*----------------------------------------------------------------------*\
790 * selected_polygon()
791 *
792 * Draws a polygon from points in x[] and y[].
793 * Point selected by v1..v4
794 \*----------------------------------------------------------------------*/
795
796 static void
selected_polygon(void (* fill)(PLINT,PLFLT *,PLFLT *),PLINT (* defined)(PLFLT,PLFLT),PLFLT * x,PLFLT * y,PLINT v1,PLINT v2,PLINT v3,PLINT v4)797 selected_polygon(void (*fill) (PLINT, PLFLT *, PLFLT *),
798 PLINT (*defined) (PLFLT, PLFLT),
799 PLFLT *x, PLFLT *y, PLINT v1, PLINT v2, PLINT v3, PLINT v4)
800 {
801 register PLINT n = 0;
802 PLFLT xx[4], yy[4];
803
804 if (fill == NULL)
805 return;
806 if (v1 >= 0) {
807 xx[n] = x[v1];
808 yy[n++] = y[v1];
809 }
810 if (v2 >= 0) {
811 xx[n] = x[v2];
812 yy[n++] = y[v2];
813 }
814 if (v3 >= 0) {
815 xx[n] = x[v3];
816 yy[n++] = y[v3];
817 }
818 if (v4 >= 0) {
819 xx[n] = x[v4];
820 yy[n++] = y[v4];
821 }
822 exfill (fill, defined, n, xx, yy);
823 }
824
825 /*----------------------------------------------------------------------*\
826 * bisect()
827 *
828 * Find boundary recursively by bisection.
829 * (x1, y1) is in the defined region, while (x2, y2) in the undefined one.
830 * The result is returned in
831 \*----------------------------------------------------------------------*/
832
833 static void
bisect(PLINT (* defined)(PLFLT,PLFLT),PLINT niter,PLFLT xx1,PLFLT yy1,PLFLT xx2,PLFLT yy2,PLFLT * xb,PLFLT * yb)834 bisect(PLINT (*defined) (PLFLT, PLFLT), PLINT niter,
835 PLFLT xx1, PLFLT yy1, PLFLT xx2, PLFLT yy2, PLFLT* xb, PLFLT* yb)
836 {
837 PLFLT xm;
838 PLFLT ym;
839
840 if (niter == 0) {
841 *xb = xx1;
842 *yb = yy1;
843 return;
844 }
845
846 xm = (xx1 + xx2) / 2;
847 ym = (yy1 + yy2) / 2;
848
849 if (defined (xm, ym))
850 bisect (defined, niter - 1, xm, ym, xx2, yy2, xb, yb);
851 else
852 bisect (defined, niter - 1, xx1, yy1, xm, ym, xb, yb);
853 }
854
855 /*----------------------------------------------------------------------*\
856 * exfill()
857 *
858 * Draws a polygon from points in x[] and y[] by taking into account
859 * eventual exclusions
860 \*----------------------------------------------------------------------*/
861
862 static void
exfill(void (* fill)(PLINT,PLFLT *,PLFLT *),PLINT (* defined)(PLFLT,PLFLT),int n,PLFLT * x,PLFLT * y)863 exfill(void (*fill) (PLINT, PLFLT *, PLFLT *),
864 PLINT (*defined) (PLFLT, PLFLT),
865 int n, PLFLT *x, PLFLT *y)
866 {
867 if (defined == NULL)
868
869 (*fill) (n, x, y);
870
871 else {
872 PLFLT xx[16];
873 PLFLT yy[16];
874 PLFLT xb, yb;
875 PLINT count = 0;
876 PLINT is_inside = defined (x[n-1], y[n-1]);
877 PLINT i;
878
879 for (i = 0; i < n; i++) {
880
881 if (defined(x[i], y[i])) {
882 if (!is_inside) {
883 if (i > 0)
884 bisect (defined, 10,
885 x[i], y[i], x[i-1], y[i-1], &xb, &yb);
886 else
887 bisect (defined, 10,
888 x[i], y[i], x[n-1], y[n-1], &xb, &yb);
889 xx[count] = xb;
890 yy[count++] = yb;
891 }
892 xx[count] = x[i];
893 yy[count++] = y[i];
894 is_inside = 1;
895 }
896 else {
897 if (is_inside) {
898 if (i > 0)
899 bisect (defined, 10,
900 x[i-1], y[i-1], x[i], y[i], &xb, &yb);
901 else
902 bisect (defined, 10,
903 x[n-1], y[n-1], x[i], y[i], &xb, &yb);
904 xx[count] = xb;
905 yy[count++] = yb;
906 is_inside = 0;
907 }
908 }
909 }
910
911 if (count)
912 (*fill) (count, xx, yy);
913 }
914 }
915
916 /*----------------------------------------------------------------------*\
917 * big_recl()
918 *
919 * find a big rectangle for shading
920 *
921 * 2 goals: minimize calls to (*fill)()
922 * keep ratio 1:3 for biggest rectangle
923 *
924 * only called by plshade()
925 *
926 * assumed that a 1 x 1 square already fits
927 *
928 * c[] = condition codes
929 * ny = c[1][0] == c[ny] (you know what I mean)
930 *
931 * returns ix, iy = length of rectangle in grid units
932 *
933 * ix < dx - 1
934 * iy < dy - 1
935 *
936 * If iy == 1 -> ix = 1 (so that cond code can be set to skip)
937 \*----------------------------------------------------------------------*/
938
939 #define RATIO 3
940 #define COND(x,y) cond_code[x*ny + y]
941
942 static void
big_recl(int * cond_code,register int ny,int dx,int dy,int * ix,int * iy)943 big_recl(int *cond_code, register int ny, int dx, int dy,
944 int *ix, int *iy)
945 {
946
947 int ok_x, ok_y, j;
948 register int i, x, y;
949 register int *cond;
950
951 /* ok_x = ok to expand in x direction */
952 /* x = current number of points in x direction */
953
954 ok_x = ok_y = 1;
955 x = y = 2;
956
957 while (ok_x || ok_y) {
958 #ifdef RATIO
959 if (RATIO * x <= y || RATIO * y <= x)
960 break;
961 #endif
962 if (ok_y) {
963 /* expand in vertical */
964 ok_y = 0;
965 if (y == dy)
966 continue;
967 cond = &COND(0, y);
968 for (i = 0; i < x; i++) {
969 if (*cond != OK)
970 break;
971 cond += ny;
972 }
973 if (i == x) {
974 /* row is ok */
975 y++;
976 ok_y = 1;
977 }
978 }
979 if (ok_x) {
980 if (y == 2)
981 break;
982 /* expand in x direction */
983 ok_x = 0;
984 if (x == dx)
985 continue;
986 cond = &COND(x, 0);
987 for (i = 0; i < y; i++) {
988 if (*cond++ != OK)
989 break;
990 }
991 if (i == y) {
992 /* column is OK */
993 x++;
994 ok_x = 1;
995 }
996 }
997 }
998
999 /* found the largest rectangle of 'ix' by 'iy' */
1000 *ix = --x;
1001 *iy = --y;
1002
1003 /* set condition code to UNDEF in interior of rectangle */
1004
1005 for (i = 1; i < x; i++) {
1006 cond = &COND(i, 1);
1007 for (j = 1; j < y; j++) {
1008 *cond++ = UNDEF;
1009 }
1010 }
1011 }
1012
1013 /*----------------------------------------------------------------------*\
1014 * draw_boundary()
1015 *
1016 * Draw boundaries of contour regions based on min_pts[], and max_pts[].
1017 \*----------------------------------------------------------------------*/
1018
1019 static void
draw_boundary(PLINT slope,PLFLT * x,PLFLT * y)1020 draw_boundary(PLINT slope, PLFLT *x, PLFLT *y)
1021 {
1022 int i;
1023
1024 if (pen_col_min != 0 && pen_wd_min != 0 && min_points != 0) {
1025 plcol0(pen_col_min);
1026 plwid(pen_wd_min);
1027 if (min_points == 4 && slope == 0) {
1028 /* swap points 1 and 3 */
1029 i = min_pts[1];
1030 min_pts[1] = min_pts[3];
1031 min_pts[3] = i;
1032 }
1033 pljoin(x[min_pts[0]], y[min_pts[0]], x[min_pts[1]], y[min_pts[1]]);
1034 if (min_points == 4) {
1035 pljoin(x[min_pts[2]], y[min_pts[2]], x[min_pts[3]],
1036 y[min_pts[3]]);
1037 }
1038 }
1039 if (pen_col_max != 0 && pen_wd_max != 0 && max_points != 0) {
1040 plcol0(pen_col_max);
1041 plwid(pen_wd_max);
1042 if (max_points == 4 && slope == 0) {
1043 /* swap points 1 and 3 */
1044 i = max_pts[1];
1045 max_pts[1] = max_pts[3];
1046 max_pts[3] = i;
1047 }
1048 pljoin(x[max_pts[0]], y[max_pts[0]], x[max_pts[1]], y[max_pts[1]]);
1049 if (max_points == 4) {
1050 pljoin(x[max_pts[2]], y[max_pts[2]], x[max_pts[3]],
1051 y[max_pts[3]]);
1052 }
1053 }
1054 }
1055
1056 /*----------------------------------------------------------------------*\
1057 *
1058 * plctest( &(x[0][0]), PLFLT level)
1059 * where x was defined as PLFLT x[4][4];
1060 *
1061 * determines if the contours associated with level have
1062 * positive slope or negative slope in the box:
1063 *
1064 * (2,3) (3,3)
1065 *
1066 * (2,2) (3,2)
1067 *
1068 * this is heuristic and can be changed by the user
1069 *
1070 * return 1 if positive slope
1071 * 0 if negative slope
1072 *
1073 * algorithmn:
1074 * 1st test:
1075 * find length of contours assuming positive and negative slopes
1076 * if the length of the negative slope contours is much bigger
1077 * than the positive slope, then the slope is positive.
1078 * (and vice versa)
1079 * (this test tries to minimize the length of contours)
1080 *
1081 * 2nd test:
1082 * if abs((top-right corner) - (bottom left corner)) >
1083 * abs((top-left corner) - (bottom right corner)) ) then
1084 * return negatiave slope.
1085 * (this test tries to keep the slope for different contour levels
1086 * the same)
1087 \*----------------------------------------------------------------------*/
1088
1089 #define X(a,b) (x[a*4+b])
1090 #define POSITIVE_SLOPE (PLINT) 1
1091 #define NEGATIVE_SLOPE (PLINT) 0
1092 #define RATIO_SQ 6.0
1093
1094 static PLINT
plctest(PLFLT * x,PLFLT level)1095 plctest(PLFLT *x, PLFLT level)
1096 {
1097 int i, j;
1098 double t[4], sorted[4], temp;
1099
1100 (void) level; /* pmr: make it used */
1101
1102 sorted[0] = t[0] = X(1,1);
1103 sorted[1] = t[1] = X(2,2);
1104 sorted[2] = t[2] = X(1,2);
1105 sorted[3] = t[3] = X(2,1);
1106
1107 for (j = 1; j < 4; j++) {
1108 temp = sorted[j];
1109 i = j - 1;
1110 while (i >= 0 && sorted[i] > temp) {
1111 sorted[i+1] = sorted[i];
1112 i--;
1113 }
1114 sorted[i+1] = temp;
1115 }
1116 /* sorted[0] == min */
1117
1118 /* find min contour */
1119 temp = int_val * ceil(sorted[0]/int_val);
1120 if (temp < sorted[1]) {
1121 /* one contour line */
1122 for (i = 0; i < 4; i++) {
1123 if (t[i] < temp) return i/2;
1124 }
1125 }
1126
1127 /* find max contour */
1128 temp = int_val * floor(sorted[3]/int_val);
1129 if (temp > sorted[2]) {
1130 /* one contour line */
1131 for (i = 0; i < 4; i++) {
1132 if (t[i] > temp) return i/2;
1133 }
1134 }
1135 /* nothing better to do - be consistant */
1136 return POSITIVE_SLOPE;
1137 }
1138
1139 /*----------------------------------------------------------------------*\
1140 * plctestez
1141 *
1142 * second routine - easier to use
1143 * fills in x[4][4] and calls plctest
1144 *
1145 * test location a[ix][iy] (lower left corner)
1146 \*----------------------------------------------------------------------*/
1147
1148 static PLINT
plctestez(PLFLT * a,PLINT nx,PLINT ny,PLINT ix,PLINT iy,PLFLT level)1149 plctestez(PLFLT *a, PLINT nx, PLINT ny, PLINT ix,
1150 PLINT iy, PLFLT level)
1151 {
1152
1153 PLFLT x[4][4];
1154 int i, j, ii, jj;
1155
1156 for (i = 0; i < 4; i++) {
1157 ii = ix + i - 1;
1158 ii = MAX(0, ii);
1159 ii = MIN(ii, nx - 1);
1160 for (j = 0; j < 4; j++) {
1161 jj = iy + j - 1;
1162 jj = MAX(0, jj);
1163 jj = MIN(jj, ny - 1);
1164 x[i][j] = a[ii * ny + jj];
1165 }
1166 }
1167 return plctest(&(x[0][0]), level);
1168 }
1169