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