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