Back to index

plt-scheme  4.2.1
plshade.c
Go to the documentation of this file.
00001 /* $Id: plshade.c,v 1.2 2005/03/17 21:39:21 eli Exp $
00002 
00003        Functions to shade regions on the basis of value.
00004        Can be used to shade contour plots or alone.
00005        Copyright 1993 Wesley Ebisuzaki 
00006 */
00007 
00008 /*----------------------------------------------------------------------*\
00009  * Call syntax for plshade():
00010  * 
00011  * void plshade(PLFLT *a, PLINT nx, PLINT ny, char *defined, 
00012  *     PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax, 
00013  *     PLFLT shade_min, PLFLT shade_max, 
00014  *     PLINT sh_color, PLINT sh_width, PLINT min_color, PLINT min_width,
00015  *     PLINT max_color, PLINT max_width, void (*fill)(), PLINT
00016  *     rectangular, ...)
00017  * 
00018  * arguments:
00019  * 
00020  *     PLFLT &(a[0][0])
00021  * 
00022  * Contains array to be plotted. The array must have been declared as
00023  * PLFLT a[nx][ny].  See following note on fortran-style arrays.
00024  * 
00025  *     PLINT nx, ny
00026  * 
00027  * Dimension of array "a".
00028  * 
00029  *     char &(defined[0][0])
00030  * 
00031  * Contains array of flags, 1 = data is valid, 0 = data is not valid.
00032  * This array determines which sections of the data is to be plotted.
00033  * This argument can be NULL if all the values are valid.  Must have been
00034  * declared as char defined[nx][ny].
00035  * 
00036  *     PLFLT xmin, xmax, ymin, ymax
00037  * 
00038  * Defines the "grid" coordinates.  The data a[0][0] has a position of
00039  * (xmin,ymin).
00040  * 
00041  *     void (*mapform)()
00042  * 
00043  * Transformation from `grid' coordinates to world coordinates.  This
00044  * pointer to a function can be NULL in which case the grid coordinates
00045  * are the same as the world coordinates.
00046  * 
00047  *     PLFLT shade_min, shade_max
00048  * 
00049  * Defines the interval to be shaded. If shade_max <= shade_min, plshade
00050  * does nothing.
00051  * 
00052  *     PLINT sh_cmap, PLFLT sh_color, PLINT sh_width
00053  * 
00054  * Defines color map, color map index, and width used by the fill pattern.
00055  * 
00056  *     PLINT min_color, min_width, max_color, max_width
00057  * 
00058  * Defines pen color, width used by the boundary of shaded region. The min
00059  * values are used for the shade_min boundary, and the max values are used
00060  * on the shade_max boundary.  Set color and width to zero for no plotted
00061  * boundaries.
00062  * 
00063  *     void (*fill)()
00064  * 
00065  * Routine used to fill the region.  Use plfill.  Future version of plplot
00066  * may have other fill routines.
00067  * 
00068  *     PLINT rectangular
00069  * 
00070  * Flag. Set to 1 if rectangles map to rectangles after (*mapform)() else
00071  * set to zero. If rectangular is set to 1, plshade tries to save time by
00072  * filling large rectangles.  This optimization fails if (*mapform)()
00073  * distorts the shape of rectangles.  For example a plot in polor
00074  * coordinates has to have rectangular set to zero.
00075  * 
00076  * Example mapform's:
00077  * 
00078  * Grid to world coordinate transformation.
00079  * This example goes from a r-theta to x-y for a polar plot.
00080  *
00081  * void mapform(PLINT n, PLFLT *x, PLFLT *y) {
00082  *     int i;
00083  *     double r, theta;
00084  *     for (i = 0; i < n; i++) {
00085  *         r = x[i];
00086  *         theta = y[i];
00087  *         x[i] = r*cos(theta);
00088  *         y[i] = r*sin(theta);    
00089  *     }
00090  * }
00091  * 
00092  * Grid was in cm, convert to world coordinates in inches.
00093  * Expands in x direction.
00094  *
00095  * void mapform(PLINT n, PLFLT *x, PLFLT *y) {
00096  *     int i;
00097  *     for (i = 0; i < n; i++) {
00098  *            x[i] = (1.0 / 2.5) * x[i];
00099  *            y[i] = (1.0 / 2.5) * y[i];
00100  *     }
00101  * }
00102  *
00103 \*----------------------------------------------------------------------*/
00104 
00105 #include "plplotP.h"
00106 #include <float.h>
00107 
00108 #define MISSING_MIN_DEF (PLFLT) 1.0
00109 #define MISSING_MAX_DEF (PLFLT) -1.0
00110 
00111 
00112 #define NEG  1
00113 #define POS  8
00114 #define OK   0
00115 #define UNDEF 64
00116 
00117 #define linear(val1, val2, level)  ((level - val1) / (val2 - val1))
00118 
00119 /* Global variables */
00120 
00121 static PLFLT sh_max, sh_min;
00122 static int min_points, max_points, n_point;
00123 static int min_pts[4], max_pts[4];
00124 static PLINT pen_col_min, pen_col_max;
00125 static PLINT pen_wd_min, pen_wd_max;
00126 static PLFLT int_val;
00127 
00128 /* Function prototypes */
00129 
00130 static void 
00131 set_cond(register int *cond, register PLFLT *a, register PLINT n);
00132 
00133 static int 
00134 find_interval(PLFLT a0, PLFLT a1, PLINT c0, PLINT c1, PLFLT *x);
00135 
00136 static void 
00137 poly(void (*fill) (PLINT, PLFLT *, PLFLT *),
00138      PLINT (*defined) (PLFLT, PLFLT),     
00139      PLFLT *x, PLFLT *y, PLINT v1, PLINT v2, PLINT v3, PLINT v4);
00140 
00141 static void 
00142 exfill(void (*fill) (PLINT, PLFLT *, PLFLT *),
00143        PLINT (*defined) (PLFLT, PLFLT),
00144        int n, PLFLT *x, PLFLT *y);
00145 
00146 static void 
00147 big_recl(int *cond_code, register int ny, int dx, int dy,
00148         int *ix, int *iy);
00149 
00150 static void 
00151 draw_boundary(PLINT slope, PLFLT *x, PLFLT *y);
00152 
00153 static PLINT 
00154 plctest(PLFLT *x, PLFLT level);
00155 
00156 static PLINT 
00157 plctestez(PLFLT *a, PLINT nx, PLINT ny, PLINT ix,
00158          PLINT iy, PLFLT level);
00159 
00160 static void
00161 plshade_int(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
00162        PLPointer f2eval_data,
00163        PLFLT (*c2eval) (PLINT, PLINT, PLPointer),
00164        PLPointer c2eval_data, 
00165        PLINT (*defined) (PLFLT, PLFLT),
00166        PLFLT missing_min, PLFLT missing_max,
00167        PLINT nx, PLINT ny, 
00168        PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
00169        PLFLT shade_min, PLFLT shade_max,
00170        PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
00171        PLINT min_color, PLINT min_width,
00172        PLINT max_color, PLINT max_width,
00173        void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
00174        void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
00175        PLPointer pltr_data);
00176 
00177 /*----------------------------------------------------------------------*\
00178  * plshades()
00179  *
00180  * Shade regions via a series of calls to plshade.
00181  * All arguments are the same as plshade except the following:
00182  * clevel is a pointer to an array of values representing 
00183  * the shade edge values, nlevel-1 is
00184  * the number of different shades, (nlevel is the number of shade edges),
00185  * fill_width is the pattern fill width, and cont_color and cont_width
00186  * are the color and width of the contour drawn at each shade edge.
00187  * (if cont_color <= 0 or cont_width <=0, no such contours are drawn).
00188 
00189 \*----------------------------------------------------------------------*/
00190 
00191 MZ_DLLEXPORT
00192 void c_plshades( PLFLT **a, PLINT nx, PLINT ny, PLINT (*defined) (PLFLT, PLFLT),
00193               PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
00194               PLFLT *clevel, PLINT nlevel, PLINT fill_width,
00195               PLINT cont_color, PLINT cont_width,
00196               void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
00197               void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
00198               PLPointer pltr_data )
00199 {
00200    PLFLT shade_min, shade_max, shade_color;
00201    PLINT i, init_color, init_width;
00202 
00203    for (i = 0; i < nlevel-1; i++) {
00204       shade_min = clevel[i];
00205       shade_max = clevel[i+1];
00206       shade_color = i / (PLFLT) (nlevel-2);
00207       /* The constants in order mean 
00208        * (1) color map1,
00209        * (0, 0, 0, 0) all edge effects will be done with plcont rather
00210        * than the normal plshade drawing which gets partially blocked
00211        * when sequential shading is done as in the present case */
00212       
00213       plshade(a, nx, ny, defined, xmin, xmax, ymin, ymax,
00214              shade_min, shade_max,
00215              1, shade_color, fill_width,
00216              0, 0, 0, 0,
00217              fill, rectangular, pltr, pltr_data);
00218    }
00219    if(cont_color > 0 && cont_width > 0) {
00220       init_color = plsc->icol0;
00221       init_width = plsc->width;
00222       plcol0(cont_color);
00223       plwid(cont_width);
00224       plcont(a, nx, ny, 1, nx, 1, ny, clevel, nlevel, pltr, pltr_data);
00225       plcol0(init_color);
00226       plwid(init_width);
00227    }
00228 }
00229    
00230 /*----------------------------------------------------------------------*\
00231  * plshade()
00232  *
00233  * Shade region.
00234  * This interface to plfshade() assumes the 2d function array is passed
00235  * via a (PLFLT **), and is column-dominant (normal C ordering).
00236 \*----------------------------------------------------------------------*/
00237 
00238 void c_plshade( PLFLT **a, PLINT nx, PLINT ny, PLINT (*defined) (PLFLT, PLFLT),
00239                  PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
00240                  PLFLT shade_min, PLFLT shade_max,
00241                  PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
00242                  PLINT min_color, PLINT min_width,
00243                  PLINT max_color, PLINT max_width,
00244                  void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
00245                  void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
00246                  PLPointer pltr_data )
00247 {
00248     PLfGrid2 grid;
00249 
00250     grid.f = a;
00251     grid.nx = nx;
00252     grid.ny = ny;
00253 
00254     plshade_int( plf2eval2, (PLPointer) &grid,
00255                  NULL, NULL,
00256 /*          plc2eval, (PLPointer) &cgrid,*/
00257                  defined, MISSING_MIN_DEF, MISSING_MAX_DEF, nx, ny, xmin, 
00258                  xmax, ymin, ymax, shade_min, shade_max,
00259                  sh_cmap, sh_color, sh_width,
00260                  min_color, min_width, max_color, max_width,
00261                  fill, rectangular, pltr, pltr_data );
00262 }
00263 
00264 /*----------------------------------------------------------------------*\
00265  * plshade1()
00266  *
00267  * Shade region.
00268  * This interface to plfshade() assumes the 2d function array is passed
00269  * via a (PLFLT *), and is column-dominant (normal C ordering).
00270 \*----------------------------------------------------------------------*/
00271 
00272 void c_plshade1( PLFLT *a, PLINT nx, PLINT ny, PLINT (*defined) (PLFLT, PLFLT),
00273                  PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
00274                  PLFLT shade_min, PLFLT shade_max,
00275                  PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
00276                  PLINT min_color, PLINT min_width,
00277                  PLINT max_color, PLINT max_width,
00278                  void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
00279                  void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
00280                  PLPointer pltr_data )
00281 {
00282     PLfGrid grid;
00283 
00284     grid.f = a;
00285     grid.nx = nx;
00286     grid.ny = ny;
00287 
00288     plshade_int( plf2eval, (PLPointer) &grid,
00289                  NULL, NULL,
00290 /*          plc2eval, (PLPointer) &cgrid,*/
00291                  defined, MISSING_MIN_DEF, MISSING_MAX_DEF, nx, ny, xmin, 
00292                  xmax, ymin, ymax, shade_min, shade_max,
00293                  sh_cmap, sh_color, sh_width,
00294                  min_color, min_width, max_color, max_width,
00295                  fill, rectangular, pltr, pltr_data );
00296 }
00297 
00298 /*----------------------------------------------------------------------*\
00299  * plfshade()
00300  *
00301  * Shade region.  
00302  * Array values are determined by the input function and the passed data.
00303 \*----------------------------------------------------------------------*/
00304 
00305 void 
00306 plfshade(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
00307         PLPointer f2eval_data,
00308         PLFLT (*c2eval) (PLINT, PLINT, PLPointer),
00309         PLPointer c2eval_data,
00310         PLINT nx, PLINT ny, 
00311         PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
00312         PLFLT shade_min, PLFLT shade_max,
00313         PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
00314         PLINT min_color, PLINT min_width,
00315         PLINT max_color, PLINT max_width,
00316         void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
00317         void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
00318         PLPointer pltr_data)
00319 {
00320     plshade_int(f2eval,  f2eval_data, c2eval, c2eval_data, 
00321         NULL, MISSING_MIN_DEF, MISSING_MAX_DEF,
00322         nx, ny, xmin, xmax, ymin, ymax,
00323         shade_min, shade_max, sh_cmap, sh_color, sh_width,
00324         min_color, min_width, max_color, max_width,
00325         fill, rectangular, pltr, pltr_data);
00326 }
00327 
00328 
00329 /*----------------------------------------------------------------------*\
00330  * plshade_int()
00331  *
00332  * Shade region -- this routine does all the work
00333  *
00334  * This routine is internal so the arguments can and will change.
00335  * To retain some compatibility between versions, you must go through
00336  * some stub routine!
00337  *
00338  * 4/95
00339  *
00340  * new: missing_min, missing_max
00341  *
00342  *     if data <= missing_max and data >= missing_min
00343  *       then the data will beconsidered to be missing
00344  *     this allows 2nd way to set undefined points (good for ftn)
00345  *     if missing feature is not used, set missing_max < missing_min
00346  *
00347  * parameters:
00348  *
00349  * f2eval, f2eval_data:  data to plot
00350  * c2eval, c2eval_data:  defined mask (not implimented)
00351  * defined: defined mask (old API - implimented)
00352  * missing_min, missing_max: yet another way to set data to undefined
00353  * nx, ny: array dimensions
00354  * xmin, xmax, ymin, ymax: grid coordinates
00355  * shade_min, shade_max: shade region with values between ...
00356  * sh_cmap, sh_color, sh_width: shading parameters, width is only for hatching
00357  * min_color, min_width: line parameters for boundary (minimum)
00358  * max_color, max_width: line parameters for boundary (maximum)
00359  *     set min_width == 0 and max_width == 0 for no contours
00360  * fill: fill function, set to NULL for no shading (contour plot)
00361  * rectangular: flag set to 1 if pltr() maps rectangles to rectangles
00362  *     this helps optimize the plotting
00363  * pltr: function to map from grid to plot coordinates
00364  *
00365  *
00366 \*----------------------------------------------------------------------*/
00367 
00368 static void
00369 plshade_int(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
00370        PLPointer f2eval_data,
00371        PLFLT (*c2eval) (PLINT, PLINT, PLPointer),
00372        PLPointer c2eval_data, 
00373        PLINT (*defined) (PLFLT, PLFLT),
00374        PLFLT missing_min, PLFLT missing_max,
00375        PLINT nx, PLINT ny, 
00376        PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
00377        PLFLT shade_min, PLFLT shade_max,
00378        PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
00379        PLINT min_color, PLINT min_width,
00380        PLINT max_color, PLINT max_width,
00381        void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
00382        void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
00383        PLPointer pltr_data)
00384 {
00385 
00386     PLINT init_width, n, slope = 0, ix, iy;
00387     int count, i, j, nxny;
00388     PLFLT *a, *a0, *a1, dx, dy;
00389     PLFLT x[8], y[8], xp[2], tx, ty;
00390     int *c, *c0, *c1;
00391 
00392     if (plsc->level < 3) {
00393        plabort("plfshade: window must be set up first");
00394        return;
00395     }
00396 
00397     if (nx <= 0 || ny <= 0) {
00398        plabort("plfshade: nx and ny must be positive");
00399        return;
00400     }
00401 
00402     if (shade_min >= shade_max) {
00403        plabort("plfshade: shade_max must exceed shade_min");
00404        return;
00405     }
00406 
00407     if (pltr == NULL || pltr_data == NULL)
00408        rectangular = 1;
00409 
00410     int_val = shade_max - shade_min;
00411     init_width = plsc->width;
00412 
00413     pen_col_min = min_color;
00414     pen_col_max = max_color;
00415 
00416     pen_wd_min = min_width;
00417     pen_wd_max = max_width;
00418 
00419     plstyl((PLINT) 0, NULL, NULL);
00420     plwid(sh_width); 
00421     if (fill != NULL) {
00422         switch (sh_cmap) {
00423         case 0:
00424             plcol0((PLINT) sh_color);
00425            break;
00426         case 1:
00427            plcol1(sh_color);
00428            break;
00429         default:
00430            plabort("plfshade: invalid color map selection");
00431            return;
00432         }
00433     }
00434     /* alloc space for value array, and initialize */
00435     /* This is only a temporary kludge */
00436     nxny = nx * ny;
00437     if ((a = (PLFLT *) malloc(nxny * sizeof(PLFLT))) == NULL) {
00438        plabort("plfshade: unable to allocate memory for value array");
00439        return;
00440     }
00441 
00442     for (ix = 0; ix < nx; ix++) 
00443        for (iy = 0; iy < ny; iy++) 
00444            a[iy + ix*ny] = f2eval(ix, iy, f2eval_data);
00445 
00446     /* alloc space for condition codes */
00447 
00448     if ((c = (int *) malloc(nxny * sizeof(int))) == NULL) {
00449        plabort("plfshade: unable to allocate memory for condition codes");
00450        free(a);
00451        return;
00452     }
00453 
00454     sh_min = shade_min;
00455     sh_max = shade_max;
00456 
00457     set_cond(c, a, nxny);
00458     dx = (xmax - xmin) / (nx - 1);
00459     dy = (ymax - ymin) / (ny - 1);
00460     a0 = a;
00461     a1 = a + ny;
00462     c0 = c;
00463     c1 = c + ny;
00464 
00465     for (ix = 0; ix < nx - 1; ix++) {
00466 
00467        for (iy = 0; iy < ny - 1; iy++) {
00468 
00469            count = c0[iy] + c0[iy + 1] + c1[iy] + c1[iy + 1];
00470 
00471            /* No filling needs to be done for these cases */
00472 
00473            if (count >= UNDEF)
00474               continue;
00475            if (count == 4 * POS)
00476               continue;
00477            if (count == 4 * NEG)
00478               continue;
00479 
00480            /* Entire rectangle can be filled */
00481 
00482            if (count == 4 * OK) {
00483               /* find biggest rectangle that fits */
00484               if (rectangular) {
00485                   big_recl(c0 + iy, ny, nx - ix, ny - iy, &i, &j);
00486               }
00487               else {
00488                   i = j = 1;
00489               }
00490               x[0] = x[1] = ix;
00491               x[2] = x[3] = ix+i;
00492               y[0] = y[3] = iy;
00493               y[1] = y[2] = iy+j;
00494 
00495               if (pltr && pltr_data) {
00496                   for (i = 0; i < 4; i++) {
00497                      (*pltr) (x[i], y[i], &tx, &ty, pltr_data);
00498                      x[i] = tx;
00499                      y[i] = ty;
00500                   }
00501               }
00502               else {
00503                   for (i = 0; i < 4; i++) {
00504                      x[i] = xmin + x[i]*dx;
00505                      y[i] = ymin + y[i]*dy;
00506                   }
00507               }
00508               if (fill != NULL)
00509                   exfill (fill, defined, (PLINT) 4, x, y);
00510               iy += j - 1;
00511               continue;
00512            }
00513 
00514            /* Only part of rectangle can be filled */
00515 
00516            n_point = min_points = max_points = 0;
00517            n = find_interval(a0[iy], a0[iy + 1], c0[iy], c0[iy + 1], xp);
00518            for (j = 0; j < n; j++) {
00519               x[j] = ix;
00520               y[j] = iy + xp[j];
00521            }
00522 
00523            i = find_interval(a0[iy + 1], a1[iy + 1],
00524                            c0[iy + 1], c1[iy + 1], xp);
00525 
00526            for (j = 0; j < i; j++) {
00527               x[j + n] = ix + xp[j];
00528               y[j + n] = iy + 1;
00529            }
00530            n += i;
00531 
00532            i = find_interval(a1[iy + 1], a1[iy], c1[iy + 1], c1[iy], xp);
00533            for (j = 0; j < i; j++) {
00534               x[n + j] = ix + 1;
00535               y[n + j] = iy + 1 - xp[j];
00536            }
00537            n += i;
00538 
00539            i = find_interval(a1[iy], a0[iy], c1[iy], c0[iy], xp);
00540            for (j = 0; j < i; j++) {
00541               x[n + j] = ix + 1 - xp[j];
00542               y[n + j] = iy;
00543            }
00544            n += i;
00545 
00546            if (pltr && pltr_data) {
00547               for (i = 0; i < n; i++) {
00548                   (*pltr) (x[i], y[i], &tx, &ty, pltr_data);
00549                   x[i] = tx;
00550                   y[i] = ty;
00551               }
00552            }
00553            else {
00554               for (i = 0; i < n; i++) {
00555                   x[i] = xmin + x[i]*dx;
00556                   y[i] = ymin + y[i]*dy;
00557               }
00558            }
00559 
00560            if (min_points == 4)
00561               slope = plctestez(a, nx, ny, ix, iy, shade_min);
00562            if (max_points == 4)
00563               slope = plctestez(a, nx, ny, ix, iy, shade_max);
00564 
00565            /* n = number of end of line segments */
00566            /* min_points = number times shade_min meets edge */
00567            /* max_points = number times shade_max meets edge */
00568 
00569            /* special cases: check number of times a contour is in a box */
00570 
00571            switch ((min_points << 3) + max_points) {
00572              case 000:
00573              case 020:
00574              case 002:
00575              case 022:
00576               if (fill != NULL && n > 0)
00577                   exfill (fill, defined, n, x, y);
00578               break;
00579              case 040:      /* 2 contour lines in box */
00580              case 004:
00581               if (n != 6)
00582                   fprintf(stderr, "plfshade err n=%d !6", (int) n);
00583               if (slope == 1 && c0[iy] == OK) {
00584                   if (fill != NULL)
00585                      exfill (fill, defined, n, x, y);
00586               }
00587               else if (slope == 1) {
00588                   poly(fill, defined, x, y, 0, 1, 2, -1);
00589                   poly(fill, defined, x, y, 3, 4, 5, -1);
00590               }
00591               else if (c0[iy + 1] == OK) {
00592                   if (fill != NULL)
00593                      exfill (fill, defined, n, x, y);
00594               }
00595               else {
00596                   poly(fill, defined, x, y, 0, 1, 5, -1);
00597                   poly(fill, defined, x, y, 2, 3, 4, -1);
00598               }
00599               break;
00600              case 044:
00601               if (n != 8)
00602                   fprintf(stderr, "plfshade err n=%d !8", (int) n);
00603               if (slope == 1) {
00604                   poly(fill, defined, x, y, 0, 1, 2, 3);
00605                   poly(fill, defined, x, y, 4, 5, 6, 7);
00606               }
00607               else {
00608                   poly(fill, defined, x, y, 0, 1, 6, 7);
00609                   poly(fill, defined, x, y, 2, 3, 4, 5);
00610               }
00611               break;
00612              case 024:
00613              case 042:
00614               /* 3 contours */
00615               if (n != 7)
00616                   fprintf(stderr, "plfshade err n=%d !7", (int) n);
00617 
00618               if ((c0[iy] == OK || c1[iy+1] == OK) && slope == 1) {
00619                   if (fill != NULL)
00620                       exfill (fill, defined, n, x, y);
00621               }
00622               else if ((c0[iy+1] == OK || c1[iy] == OK) && slope == 0) {
00623                   if (fill !=NULL)
00624                       exfill (fill, defined, n, x, y);
00625               }
00626 
00627               else if (c0[iy] == OK) {
00628                   poly(fill, defined, x, y, 0, 1, 6, -1);
00629                   poly(fill, defined, x, y, 2, 3, 4, 5);
00630               }
00631               else if (c0[iy+1] == OK) {
00632                   poly(fill, defined, x, y, 0, 1, 2, -1);
00633                   poly(fill, defined, x, y, 3, 4, 5, 6);
00634               }
00635               else if (c1[iy+1] == OK) {
00636                   poly(fill, defined, x, y, 0, 1, 5, 6);
00637                   poly(fill, defined, x, y, 2, 3, 4, -1);
00638               }
00639               else if (c1[iy] == OK) {
00640                   poly(fill, defined, x, y, 0, 1, 2, 3);
00641                   poly(fill, defined, x, y, 4, 5, 6, -1);
00642               }
00643               else {
00644                   fprintf(stderr, "plfshade err logic case 024:042\n");
00645               }
00646               break;
00647              default:
00648               fprintf(stderr, "prog err switch\n");
00649               break;
00650            }
00651            draw_boundary(slope, x, y);
00652 
00653            if (fill != NULL) {
00654                plwid(sh_width);
00655               if (sh_cmap == 0) plcol0((PLINT) sh_color);
00656               else if (sh_cmap == 1) plcol1(sh_color);
00657            }
00658        }
00659 
00660        a0 = a1;
00661        c0 = c1;
00662        a1 += ny;
00663        c1 += ny;
00664     }
00665 
00666     free(c);
00667     free(a);
00668     plwid(init_width);
00669 }
00670 
00671 /*----------------------------------------------------------------------*\
00672  * set_cond()
00673  *
00674  * Fills out condition code array.
00675 \*----------------------------------------------------------------------*/
00676 
00677 static void 
00678 set_cond(register int *cond, register PLFLT *a, register PLINT n)
00679 {
00680     while (n--) {
00681         if (*a < sh_min)
00682            *cond++ = NEG;
00683        else if (*a > sh_max)
00684            *cond++ = POS;
00685        else
00686            *cond++ = OK;
00687        a++;
00688     }
00689 }
00690 
00691 /*----------------------------------------------------------------------*\
00692  * find_interval()
00693  *
00694  * Two points x(0) = a0, (condition code c0) x(1) = a1, (condition code c1)
00695  * return interval on the line that are shaded
00696  * 
00697  * returns 0 : no points to be shaded 1 : x[0] <= x < 1 is the interval 2 :
00698  * x[0] <= x <= x[1] < 1 interval to be shaded n_point, max_points,
00699  * min_points are incremented location of min/max_points are stored 
00700 \*----------------------------------------------------------------------*/
00701 
00702 static int 
00703 find_interval(PLFLT a0, PLFLT a1, PLINT c0, PLINT c1, PLFLT *x)
00704 {
00705     register int n;
00706 
00707     n = 0;
00708     if (c0 == OK) {
00709        x[n++] = 0.0;
00710        n_point++;
00711     }
00712     if (c0 == c1)
00713        return n;
00714 
00715     if (c0 == NEG || c1 == POS) {
00716        if (c0 == NEG) {
00717            x[n++] = linear(a0, a1, sh_min);
00718            min_pts[min_points++] = n_point++;
00719        }
00720        if (c1 == POS) {
00721            x[n++] = linear(a0, a1, sh_max);
00722            max_pts[max_points++] = n_point++;
00723        }
00724     }
00725     if (c0 == POS || c1 == NEG) {
00726        if (c0 == POS) {
00727            x[n++] = linear(a0, a1, sh_max);
00728            max_pts[max_points++] = n_point++;
00729        }
00730        if (c1 == NEG) {
00731            x[n++] = linear(a0, a1, sh_min);
00732            min_pts[min_points++] = n_point++;
00733        }
00734     }
00735     return n;
00736 }
00737 
00738 /*----------------------------------------------------------------------*\
00739  * poly()
00740  *
00741  * Draws a polygon from points in x[] and y[].
00742  * Point selected by v1..v4 
00743 \*----------------------------------------------------------------------*/
00744 
00745 static void 
00746 poly(void (*fill) (PLINT, PLFLT *, PLFLT *),
00747      PLINT (*defined) (PLFLT, PLFLT),
00748      PLFLT *x, PLFLT *y, PLINT v1, PLINT v2, PLINT v3, PLINT v4)
00749 {
00750     register PLINT n = 0;
00751     PLFLT xx[4], yy[4];
00752 
00753     if (fill == NULL)
00754        return;
00755     if (v1 >= 0) {
00756        xx[n] = x[v1];
00757        yy[n++] = y[v1];
00758     }
00759     if (v2 >= 0) {
00760        xx[n] = x[v2];
00761        yy[n++] = y[v2];
00762     }
00763     if (v3 >= 0) {
00764        xx[n] = x[v3];
00765        yy[n++] = y[v3];
00766     }
00767     if (v4 >= 0) {
00768        xx[n] = x[v4];
00769        yy[n++] = y[v4];
00770     }
00771     exfill (fill, defined, n, xx, yy);
00772 }
00773 
00774 /*----------------------------------------------------------------------*\
00775  * bisect()
00776  *
00777  * Find boundary recursively by bisection.  
00778  * (x1, y1) is in the defined region, while (x2, y2) in the undefined one.
00779  * The result is returned in 
00780 \*----------------------------------------------------------------------*/
00781 
00782 static void 
00783 bisect(PLINT (*defined) (PLFLT, PLFLT), PLINT niter,
00784        PLFLT x1, PLFLT y1, PLFLT x2, PLFLT y2, PLFLT* xb, PLFLT* yb)
00785 {
00786     PLFLT xm;
00787     PLFLT ym;
00788 
00789     if (niter == 0) {
00790         *xb = x1;
00791         *yb = y1;
00792         return;
00793     }
00794 
00795     xm = (x1 + x2) / 2;
00796     ym = (y1 + y2) / 2;
00797 
00798     if (defined (xm, ym))
00799       bisect (defined, niter - 1, xm, ym, x2, y2, xb, yb);
00800     else
00801       bisect (defined, niter - 1, x1, y1, xm, ym, xb, yb);      
00802 }
00803 
00804 /*----------------------------------------------------------------------*\
00805  * exfill()
00806  *
00807  * Draws a polygon from points in x[] and y[] by taking into account 
00808  * eventual exclusions
00809 \*----------------------------------------------------------------------*/
00810 
00811 static void 
00812 exfill(void (*fill) (PLINT, PLFLT *, PLFLT *),
00813        PLINT (*defined) (PLFLT, PLFLT),
00814        int n, PLFLT *x, PLFLT *y)
00815 {
00816     if (defined == NULL)
00817 
00818         (*fill) (n, x, y);
00819 
00820     else {
00821         PLFLT xx[16];  
00822        PLFLT yy[16];
00823        PLFLT xb, yb;  
00824        PLINT count = 0;
00825        PLINT is_inside = defined (x[n-1], y[n-1]);
00826        PLINT i;
00827 
00828         for (i = 0; i < n; i++) {
00829 
00830            if (defined(x[i], y[i])) {
00831                if (!is_inside) {
00832                   if (i > 0)
00833                       bisect (defined, 10,
00834                             x[i], y[i], x[i-1], y[i-1], &xb, &yb);
00835                   else
00836                       bisect (defined, 10,
00837                             x[i], y[i], x[n-1], y[n-1], &xb, &yb);
00838                   xx[count] = xb;
00839                   yy[count++] = yb;
00840               }
00841               xx[count] = x[i];
00842               yy[count++] = y[i];
00843               is_inside = 1;
00844            }
00845            else {
00846                if (is_inside) {
00847                   if (i > 0)
00848                       bisect (defined, 10,
00849                             x[i-1], y[i-1], x[i], y[i], &xb, &yb);
00850                   else
00851                       bisect (defined, 10,
00852                             x[n-1], y[n-1], x[i], y[i], &xb, &yb);
00853                   xx[count] = xb;
00854                   yy[count++] = yb;
00855                   is_inside = 0;
00856               }
00857            }
00858        }
00859        
00860        if (count)
00861            (*fill) (count, xx, yy);
00862     }
00863 }
00864 
00865 /*----------------------------------------------------------------------*\
00866  * big_recl()
00867  *
00868  * find a big rectangle for shading
00869  *
00870  * 2 goals: minimize calls to (*fill)()
00871  *  keep ratio 1:3 for biggest rectangle
00872  *
00873  * only called by plshade()
00874  *
00875  * assumed that a 1 x 1 square already fits
00876  *
00877  * c[] = condition codes
00878  * ny = c[1][0] == c[ny]  (you know what I mean)
00879  *
00880  * returns ix, iy = length of rectangle in grid units
00881  *
00882  * ix < dx - 1
00883  * iy < dy - 1
00884  *
00885  * If iy == 1 -> ix = 1 (so that cond code can be set to skip)
00886 \*----------------------------------------------------------------------*/
00887 
00888 #define RATIO 3
00889 #define COND(x,y) cond_code[x*ny + y]
00890 
00891 static void 
00892 big_recl(int *cond_code, register int ny, int dx, int dy,
00893         int *ix, int *iy)
00894 {
00895 
00896     int ok_x, ok_y, j;
00897     register int i, x, y;
00898     register int *cond;
00899 
00900     /* ok_x = ok to expand in x direction */
00901     /* x = current number of points in x direction */
00902 
00903     ok_x = ok_y = 1;
00904     x = y = 2;
00905 
00906     while (ok_x || ok_y) {
00907 #ifdef RATIO
00908        if (RATIO * x <= y || RATIO * y <= x)
00909            break;
00910 #endif
00911        if (ok_y) {
00912            /* expand in vertical */
00913            ok_y = 0;
00914            if (y == dy)
00915               continue;
00916            cond = &COND(0, y);
00917            for (i = 0; i < x; i++) {
00918               if (*cond != OK)
00919                   break;
00920               cond += ny;
00921            }
00922            if (i == x) {
00923               /* row is ok */
00924               y++;
00925               ok_y = 1;
00926            }
00927        }
00928        if (ok_x) {
00929            if (y == 2)
00930               break;
00931            /* expand in x direction */
00932            ok_x = 0;
00933            if (x == dx)
00934               continue;
00935            cond = &COND(x, 0);
00936            for (i = 0; i < y; i++) {
00937               if (*cond++ != OK)
00938                   break;
00939            }
00940            if (i == y) {
00941               /* column is OK */
00942               x++;
00943               ok_x = 1;
00944            }
00945        }
00946     }
00947 
00948     /* found the largest rectangle of 'ix' by 'iy' */
00949     *ix = --x;
00950     *iy = --y;
00951 
00952     /* set condition code to UNDEF in interior of rectangle */
00953 
00954     for (i = 1; i < x; i++) {
00955        cond = &COND(i, 1);
00956        for (j = 1; j < y; j++) {
00957            *cond++ = UNDEF;
00958        }
00959     }
00960 }
00961 
00962 /*----------------------------------------------------------------------*\
00963  * draw_boundary()
00964  *
00965  * Draw boundaries of contour regions based on min_pts[], and max_pts[].
00966 \*----------------------------------------------------------------------*/
00967 
00968 static void 
00969 draw_boundary(PLINT slope, PLFLT *x, PLFLT *y)
00970 {
00971     int i;
00972 
00973     if (pen_col_min != 0 && pen_wd_min != 0 && min_points != 0) {
00974        plcol0(pen_col_min);
00975        plwid(pen_wd_min);
00976        if (min_points == 4 && slope == 0) {
00977            /* swap points 1 and 3 */
00978            i = min_pts[1];
00979            min_pts[1] = min_pts[3];
00980            min_pts[3] = i;
00981        }
00982        pljoin(x[min_pts[0]], y[min_pts[0]], x[min_pts[1]], y[min_pts[1]]);
00983        if (min_points == 4) {
00984            pljoin(x[min_pts[2]], y[min_pts[2]], x[min_pts[3]],
00985                  y[min_pts[3]]);
00986        }
00987     }
00988     if (pen_col_max != 0 && pen_wd_max != 0 && max_points != 0) {
00989        plcol0(pen_col_max);
00990        plwid(pen_wd_max);
00991        if (max_points == 4 && slope == 0) {
00992            /* swap points 1 and 3 */
00993            i = max_pts[1];
00994            max_pts[1] = max_pts[3];
00995            max_pts[3] = i;
00996        }
00997        pljoin(x[max_pts[0]], y[max_pts[0]], x[max_pts[1]], y[max_pts[1]]);
00998        if (max_points == 4) {
00999            pljoin(x[max_pts[2]], y[max_pts[2]], x[max_pts[3]],
01000                  y[max_pts[3]]);
01001        }
01002     }
01003 }
01004 
01005 /*----------------------------------------------------------------------*\
01006  *
01007  * plctest( &(x[0][0]), PLFLT level)
01008  * where x was defined as PLFLT x[4][4];
01009  *
01010  * determines if the contours associated with level have
01011  * positive slope or negative slope in the box:
01012  *
01013  *  (2,3)   (3,3)
01014  *
01015  *  (2,2)   (3,2)
01016  *
01017  * this is heuristic and can be changed by the user
01018  *
01019  * return 1 if positive slope
01020  *        0 if negative slope
01021  *
01022  * algorithmn:
01023  *      1st test:
01024  *     find length of contours assuming positive and negative slopes
01025  *      if the length of the negative slope contours is much bigger
01026  *     than the positive slope, then the slope is positive.
01027  *      (and vice versa)
01028  *      (this test tries to minimize the length of contours)
01029  *
01030  *      2nd test:
01031  *      if abs((top-right corner) - (bottom left corner)) >
01032  *        abs((top-left corner) - (bottom right corner)) ) then
01033  *            return negatiave slope.
01034  *      (this test tries to keep the slope for different contour levels
01035  *     the same)
01036 \*----------------------------------------------------------------------*/
01037 
01038 #define X(a,b) (x[a*4+b])
01039 #define POSITIVE_SLOPE (PLINT) 1
01040 #define NEGATIVE_SLOPE (PLINT) 0
01041 #define RATIO_SQ 6.0
01042 
01043 static PLINT 
01044 plctest(PLFLT *x, PLFLT level)
01045 {
01046     int i, j;
01047     double t[4], sorted[4], temp;
01048 
01049     sorted[0] = t[0] = X(1,1);
01050     sorted[1] = t[1] = X(2,2);
01051     sorted[2] = t[2] = X(1,2);
01052     sorted[3] = t[3] = X(2,1);
01053 
01054     for (j = 1; j < 4; j++) {
01055        temp = sorted[j];
01056        i = j - 1;
01057        while (i >= 0 && sorted[i] > temp) {
01058            sorted[i+1] = sorted[i];
01059            i--;
01060        }
01061        sorted[i+1] = temp;
01062     }
01063     /* sorted[0] == min */
01064 
01065     /* find min contour */
01066     temp = int_val * ceil(sorted[0]/int_val);
01067     if (temp < sorted[1]) {
01068        /* one contour line */
01069        for (i = 0; i < 4; i++) {
01070            if (t[i] < temp) return i/2;
01071        }
01072     }
01073        
01074     /* find max contour */
01075     temp = int_val * floor(sorted[3]/int_val);
01076     if (temp > sorted[2]) {
01077        /* one contour line */
01078        for (i = 0; i < 4; i++) {
01079            if (t[i] > temp) return i/2;
01080        }
01081     }
01082     /* nothing better to do - be consistant */
01083     return POSITIVE_SLOPE;
01084 }
01085 
01086 /*----------------------------------------------------------------------*\
01087  * plctestez
01088  *
01089  * second routine - easier to use
01090  * fills in x[4][4] and calls plctest
01091  *
01092  * test location a[ix][iy] (lower left corner)
01093 \*----------------------------------------------------------------------*/
01094 
01095 static PLINT 
01096 plctestez(PLFLT *a, PLINT nx, PLINT ny, PLINT ix,
01097          PLINT iy, PLFLT level)
01098 {
01099 
01100     PLFLT x[4][4];
01101     int i, j, ii, jj;
01102 
01103     for (i = 0; i < 4; i++) {
01104        ii = ix + i - 1;
01105        ii = MAX(0, ii);
01106        ii = MIN(ii, nx - 1);
01107        for (j = 0; j < 4; j++) {
01108            jj = iy + j - 1;
01109            jj = MAX(0, jj);
01110            jj = MIN(jj, ny - 1);
01111            x[i][j] = a[ii * ny + jj];
01112        }
01113     }
01114     return plctest(&(x[0][0]), level);
01115 }