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