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