PLplot
5.10.0
|
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 }