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