PLplot
5.10.0
|
00001 // Contour plotter. 00002 // 00003 // Copyright (C) 1995, 2000, 2001 Maurice LeBrun 00004 // Copyright (C) 2000, 2002 Joao Cardoso 00005 // Copyright (C) 2000-2014 Alan W. Irwin 00006 // Copyright (C) 2004 Andrew Ross 00007 // 00008 // This file is part of PLplot. 00009 // 00010 // PLplot is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Library General Public License as published 00012 // by the Free Software Foundation; either version 2 of the License, or 00013 // (at your option) any later version. 00014 // 00015 // PLplot is distributed in the hope that it will be useful, 00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00018 // GNU Library General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Library General Public License 00021 // along with PLplot; if not, write to the Free Software 00022 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00023 // 00024 00025 #include "plplotP.h" 00026 00027 #ifdef MSDOS 00028 #pragma optimize("",off) 00029 #endif 00030 00031 // Static function prototypes. 00032 00033 static void 00034 plcntr( PLFLT ( *plf2eval )( PLINT, PLINT, PLPointer ), 00035 PLPointer plf2eval_data, 00036 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00037 PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts, 00038 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00039 PLPointer pltr_data ); 00040 00041 static void 00042 pldrawcn( PLFLT ( *plf2eval )( PLINT, PLINT, PLPointer ), 00043 PLPointer plf2eval_data, 00044 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00045 PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow, 00046 PLFLT lastx, PLFLT lasty, PLINT startedge, 00047 PLINT **ipts, PLFLT *distance, PLINT *lastindex, 00048 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00049 PLPointer pltr_data ); 00050 00051 static void 00052 plfloatlabel( PLFLT value, char *string, PLINT len ); 00053 00054 static PLFLT 00055 plP_pcwcx( PLINT x ); 00056 00057 static PLFLT 00058 plP_pcwcy( PLINT y ); 00059 00060 static void 00061 pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex ); 00062 00063 // Error flag for aborts 00064 00065 static int error; 00066 00067 //************************************** 00068 // 00069 // Defaults for contour label printing. 00070 // 00071 //************************************** 00072 00073 // Font height for contour labels (normalized) 00074 static PLFLT 00075 contlabel_size = 0.3; 00076 00077 // Offset of label from contour line (if set to 0.0, labels are printed on the lines). 00078 static PLFLT 00079 contlabel_offset = 0.006; 00080 00081 // Spacing parameter for contour labels 00082 static PLFLT 00083 contlabel_space = 0.1; 00084 00085 // Activate labels, default off 00086 static PLINT 00087 contlabel_active = 0; 00088 00089 // If the contour label exceed 10^(limexp) or 10^(-limexp), the exponential format is used 00090 static PLINT 00091 limexp = 4; 00092 00093 // Number of significant digits 00094 static PLINT 00095 sigprec = 2; 00096 00097 //******* contour lines storage *************************** 00098 00099 static CONT_LEVEL *startlev = NULL; 00100 static CONT_LEVEL *currlev; 00101 static CONT_LINE *currline; 00102 00103 static int cont3d = 0; 00104 00105 static CONT_LINE * 00106 alloc_line( void ) 00107 { 00108 CONT_LINE *line; 00109 00110 if ( ( line = (CONT_LINE *) malloc( sizeof ( CONT_LINE ) ) ) == NULL ) 00111 { 00112 plexit( "alloc_line: Insufficient memory" ); 00113 } 00114 00115 line->x = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) ); 00116 line->y = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) ); 00117 00118 if ( ( line->x == NULL ) || ( line->y == NULL ) ) 00119 { 00120 plexit( "alloc_line: Insufficient memory" ); 00121 } 00122 00123 line->npts = 0; 00124 line->next = NULL; 00125 00126 return line; 00127 } 00128 00129 static CONT_LEVEL * 00130 alloc_level( PLFLT level ) 00131 { 00132 CONT_LEVEL *node; 00133 00134 if ( ( node = (CONT_LEVEL *) malloc( sizeof ( CONT_LEVEL ) ) ) == NULL ) 00135 { 00136 plexit( "alloc_level: Insufficient memory" ); 00137 } 00138 node->level = level; 00139 node->next = NULL; 00140 node->line = alloc_line( ); 00141 00142 return node; 00143 } 00144 00145 static void 00146 realloc_line( CONT_LINE *line ) 00147 { 00148 if ( ( ( line->x = (PLFLT *) realloc( line->x, 00149 (size_t) ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) || 00150 ( ( line->y = (PLFLT *) realloc( line->y, 00151 (size_t) ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) ) 00152 plexit( "realloc_line: Insufficient memory" ); 00153 } 00154 00155 00156 // new contour level 00157 static void 00158 cont_new_store( PLFLT level ) 00159 { 00160 if ( cont3d ) 00161 { 00162 if ( startlev == NULL ) 00163 { 00164 startlev = alloc_level( level ); 00165 currlev = startlev; 00166 } 00167 else 00168 { 00169 currlev->next = alloc_level( level ); 00170 currlev = currlev->next; 00171 } 00172 currline = currlev->line; 00173 } 00174 } 00175 00176 void 00177 cont_clean_store( CONT_LEVEL *ct ) 00178 { 00179 CONT_LINE *tline, *cline; 00180 CONT_LEVEL *tlev, *clevel; 00181 00182 if ( ct != NULL ) 00183 { 00184 clevel = ct; 00185 00186 do 00187 { 00188 cline = clevel->line; 00189 do 00190 { 00191 #ifdef CONT_PLOT_DEBUG // for 2D plots. For 3D plots look at plot3.c:plotsh3di() 00192 plP_movwor( cline->x[0], cline->y[0] ); 00193 for ( j = 1; j < cline->npts; j++ ) 00194 plP_drawor( cline->x[j], cline->y[j] ); 00195 #endif 00196 tline = cline->next; 00197 free( cline->x ); 00198 free( cline->y ); 00199 free( cline ); 00200 cline = tline; 00201 } 00202 while ( cline != NULL ); 00203 tlev = clevel->next; 00204 free( clevel ); 00205 clevel = tlev; 00206 } 00207 while ( clevel != NULL ); 00208 startlev = NULL; 00209 } 00210 } 00211 00212 static void 00213 cont_xy_store( PLFLT xx, PLFLT yy ) 00214 { 00215 if ( cont3d ) 00216 { 00217 PLINT pts = currline->npts; 00218 00219 if ( pts % LINE_ITEMS == 0 ) 00220 realloc_line( currline ); 00221 00222 currline->x[pts] = xx; 00223 currline->y[pts] = yy; 00224 currline->npts++; 00225 } 00226 else 00227 plP_drawor( xx, yy ); 00228 } 00229 00230 static void 00231 cont_mv_store( PLFLT xx, PLFLT yy ) 00232 { 00233 if ( cont3d ) 00234 { 00235 if ( currline->npts != 0 ) // not an empty list, allocate new 00236 { 00237 currline->next = alloc_line( ); 00238 currline = currline->next; 00239 } 00240 00241 // and fill first element 00242 currline->x[0] = xx; 00243 currline->y[0] = yy; 00244 currline->npts = 1; 00245 } 00246 else 00247 plP_movwor( xx, yy ); 00248 } 00249 00250 // small routine to set offset and spacing of contour labels, see desciption above 00251 void c_pl_setcontlabelparam( PLFLT offset, PLFLT size, PLFLT spacing, PLINT active ) 00252 { 00253 contlabel_offset = offset; 00254 contlabel_size = size; 00255 contlabel_space = spacing; 00256 contlabel_active = active; 00257 } 00258 00259 // small routine to set the format of the contour labels, description of limexp and prec see above 00260 void c_pl_setcontlabelformat( PLINT lexp, PLINT sigdig ) 00261 { 00262 limexp = lexp; 00263 sigprec = sigdig; 00264 } 00265 00266 static void pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex ) 00267 { 00268 PLFLT delta_x, delta_y; 00269 PLINT currx_old, curry_old; 00270 00271 delta_x = plP_pcdcx( plsc->currx ) - plP_pcdcx( plP_wcpcx( tpx ) ); 00272 delta_y = plP_pcdcy( plsc->curry ) - plP_pcdcy( plP_wcpcy( tpy ) ); 00273 00274 currx_old = plsc->currx; 00275 curry_old = plsc->curry; 00276 00277 *distance += sqrt( delta_x * delta_x + delta_y * delta_y ); 00278 00279 plP_drawor( tpx, tpy ); 00280 00281 if ( (int) ( fabs( *distance / contlabel_space ) ) > *lastindex ) 00282 { 00283 PLFLT scale, vec_x, vec_y, mx, my, dev_x, dev_y, off_x, off_y; 00284 00285 vec_x = tpx - plP_pcwcx( currx_old ); 00286 vec_y = tpy - plP_pcwcy( curry_old ); 00287 00288 // Ensure labels appear the right way up 00289 if ( vec_x < 0 ) 00290 { 00291 vec_x = -vec_x; 00292 vec_y = -vec_y; 00293 } 00294 00295 mx = (double) plsc->wpxscl / (double) plsc->phyxlen; 00296 my = (double) plsc->wpyscl / (double) plsc->phyylen; 00297 00298 dev_x = -my * vec_y / mx; 00299 dev_y = mx * vec_x / my; 00300 00301 scale = sqrt( ( mx * mx * dev_x * dev_x + my * my * dev_y * dev_y ) / 00302 ( contlabel_offset * contlabel_offset ) ); 00303 00304 off_x = dev_x / scale; 00305 off_y = dev_y / scale; 00306 00307 plptex( tpx + off_x, tpy + off_y, vec_x, vec_y, 0.5, flabel ); 00308 plP_movwor( tpx, tpy ); 00309 ( *lastindex )++; 00310 } 00311 else 00312 plP_movwor( tpx, tpy ); 00313 } 00314 00315 00316 // Format contour labels. Arguments: 00317 // value: floating point number to be formatted 00318 // string: the formatted label, plptex must be called with it to actually 00319 // print the label 00320 // 00321 00322 static void plfloatlabel( PLFLT value, char *string, PLINT len ) 00323 { 00324 PLINT setpre, precis; 00325 // form[10] gives enough space for all non-malicious formats. 00326 // tmpstring[15] gives enough room for 3 digits in a negative exponent 00327 // or 4 digits in a positive exponent + null termination. That 00328 // should be enough for all non-malicious use. 00329 // Obviously there are security issues here that 00330 // should be addressed as well. 00331 // 00332 #define FORM_LEN 10 00333 #define TMPSTRING_LEN 15 00334 char form[FORM_LEN], tmpstring[TMPSTRING_LEN]; 00335 PLINT exponent = 0; 00336 PLFLT mant, tmp; 00337 00338 PLINT prec = sigprec; 00339 00340 plP_gprec( &setpre, &precis ); 00341 00342 if ( setpre ) 00343 prec = precis; 00344 00345 if ( value > 0.0 ) 00346 tmp = log10( value ); 00347 else if ( value < 0.0 ) 00348 tmp = log10( -value ); 00349 else 00350 tmp = 0; 00351 00352 if ( tmp >= 0.0 ) 00353 exponent = (int) tmp; 00354 else if ( tmp < 0.0 ) 00355 { 00356 tmp = -tmp; 00357 if ( floor( tmp ) < tmp ) 00358 exponent = -(int) ( floor( tmp ) + 1.0 ); 00359 else 00360 exponent = -(int) ( floor( tmp ) ); 00361 } 00362 00363 mant = value / pow( 10.0, exponent ); 00364 00365 if ( mant != 0.0 ) 00366 mant = (int) ( mant * pow( 10.0, prec - 1 ) + 0.5 * mant / fabs( mant ) ) / pow( 10.0, prec - 1 ); 00367 00368 snprintf( form, FORM_LEN, "%%.%df", prec - 1 ); 00369 snprintf( string, (size_t) len, form, mant ); 00370 snprintf( tmpstring, TMPSTRING_LEN, "#(229)10#u%d", exponent ); 00371 strncat( string, tmpstring, (size_t) len - strlen( string ) - 1 ); 00372 00373 if ( abs( exponent ) < limexp || value == 0.0 ) 00374 { 00375 value = pow( 10.0, exponent ) * mant; 00376 00377 if ( exponent >= 0 ) 00378 prec = prec - 1 - exponent; 00379 else 00380 prec = prec - 1 + abs( exponent ); 00381 00382 if ( prec < 0 ) 00383 prec = 0; 00384 00385 snprintf( form, FORM_LEN, "%%.%df", (int) prec ); 00386 snprintf( string, (size_t) len, form, value ); 00387 } 00388 } 00389 00390 // physical coords (x) to world coords 00391 00392 static PLFLT 00393 plP_pcwcx( PLINT x ) 00394 { 00395 return ( ( x - plsc->wpxoff ) / plsc->wpxscl ); 00396 } 00397 00398 // physical coords (y) to world coords 00399 00400 static PLFLT 00401 plP_pcwcy( PLINT y ) 00402 { 00403 return ( ( y - plsc->wpyoff ) / plsc->wpyscl ); 00404 } 00405 00406 //-------------------------------------------------------------------------- 00407 // plf2eval1() 00408 // 00409 // Does a lookup from a 2d function array. Array is of type (PLFLT **), 00410 // and is column dominant (normal C ordering). 00411 //-------------------------------------------------------------------------- 00412 00413 PLFLT 00414 plf2eval1( PLINT ix, PLINT iy, PLPointer plf2eval_data ) 00415 { 00416 PLFLT value; 00417 PLFLT **z = (PLFLT **) plf2eval_data; 00418 00419 value = z[ix][iy]; 00420 00421 return value; 00422 } 00423 00424 //-------------------------------------------------------------------------- 00425 // plf2eval2() 00426 // 00427 // Does a lookup from a 2d function array. plf2eval_data is treated as type 00428 // (PLfGrid2 *). 00429 //-------------------------------------------------------------------------- 00430 00431 PLFLT 00432 plf2eval2( PLINT ix, PLINT iy, PLPointer plf2eval_data ) 00433 { 00434 PLFLT value; 00435 PLfGrid2 *grid = (PLfGrid2 *) plf2eval_data; 00436 00437 value = grid->f[ix][iy]; 00438 00439 return value; 00440 } 00441 00442 //-------------------------------------------------------------------------- 00443 // plf2eval() 00444 // 00445 // Does a lookup from a 2d function array. Array is of type (PLFLT *), and 00446 // is column dominant (normal C ordering). You MUST fill the ny maximum 00447 // array index entry in the PLfGrid struct. 00448 //-------------------------------------------------------------------------- 00449 00450 PLFLT 00451 plf2eval( PLINT ix, PLINT iy, PLPointer plf2eval_data ) 00452 { 00453 PLFLT value; 00454 PLfGrid *grid = (PLfGrid *) plf2eval_data; 00455 00456 value = grid->f[ix * grid->ny + iy]; 00457 00458 return value; 00459 } 00460 00461 //-------------------------------------------------------------------------- 00462 // plf2evalr() 00463 // 00464 // Does a lookup from a 2d function array. Array is of type (PLFLT *), and 00465 // is row dominant (Fortran ordering). You MUST fill the nx maximum array 00466 // index entry in the PLfGrid struct. 00467 //-------------------------------------------------------------------------- 00468 00469 PLFLT 00470 plf2evalr( PLINT ix, PLINT iy, PLPointer plf2eval_data ) 00471 { 00472 PLFLT value; 00473 PLfGrid *grid = (PLfGrid *) plf2eval_data; 00474 00475 value = grid->f[ix + iy * grid->nx]; 00476 00477 return value; 00478 } 00479 00480 //-------------------------------------------------------------------------- 00481 // 00482 // cont_store: 00483 // 00484 // Draw contour lines in memory. 00485 // cont_clean_store() must be called after use to release allocated memory. 00486 // 00487 //-------------------------------------------------------------------------- 00488 00489 void 00490 cont_store( const PLFLT * const *f, PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00491 PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel, 00492 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00493 PLPointer pltr_data, 00494 CONT_LEVEL **contour ) 00495 { 00496 cont3d = 1; 00497 00498 plcont( f, nx, ny, kx, lx, ky, ly, clevel, nlevel, 00499 pltr, pltr_data ); 00500 00501 *contour = startlev; 00502 cont3d = 0; 00503 } 00504 00505 //-------------------------------------------------------------------------- 00506 // void plcont() 00507 // 00508 // Draws a contour plot from data in f(nx,ny). Is just a front-end to 00509 // plfcont, with a particular choice for f2eval and f2eval_data. 00510 //-------------------------------------------------------------------------- 00511 00512 void 00513 c_plcont( const PLFLT * const *f, PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00514 PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel, 00515 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00516 PLPointer pltr_data ) 00517 { 00518 if ( pltr == NULL ) 00519 { 00520 // If pltr is undefined, abort with an error. 00521 plabort( "plcont: The pltr callback must be defined" ); 00522 return; 00523 } 00524 00525 plfcont( plf2eval1, (PLPointer) f, 00526 nx, ny, kx, lx, ky, ly, clevel, nlevel, 00527 pltr, pltr_data ); 00528 } 00529 00530 //-------------------------------------------------------------------------- 00531 // void plfcont() 00532 // 00533 // Draws a contour plot using the function evaluator f2eval and data stored 00534 // by way of the f2eval_data pointer. This allows arbitrary organizations 00535 // of 2d array data to be used. 00536 // 00537 // The subrange of indices used for contouring is kx to lx in the x 00538 // direction and from ky to ly in the y direction. The array of contour 00539 // levels is clevel(nlevel), and "pltr" is the name of a function which 00540 // transforms array indicies into world coordinates. 00541 // 00542 // Note that the fortran-like minimum and maximum indices (kx, lx, ky, ly) 00543 // are translated into more C-like ones. I've only kept them as they are 00544 // for the plcontf() argument list because of backward compatibility. 00545 //-------------------------------------------------------------------------- 00546 00547 void 00548 plfcont( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ), 00549 PLPointer f2eval_data, 00550 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00551 PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel, 00552 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00553 PLPointer pltr_data ) 00554 { 00555 PLINT i, **ipts; 00556 00557 if ( kx < 1 || kx >= lx ) 00558 { 00559 plabort( "plfcont: indices must satisfy 1 <= kx <= lx <= nx" ); 00560 return; 00561 } 00562 if ( ky < 1 || ky >= ly ) 00563 { 00564 plabort( "plfcont: indices must satisfy 1 <= ky <= ly <= ny" ); 00565 return; 00566 } 00567 00568 if ( ( ipts = (PLINT **) malloc( (size_t) nx * sizeof ( PLINT * ) ) ) == NULL ) 00569 { 00570 plexit( "plfcont: Insufficient memory" ); 00571 } 00572 00573 for ( i = 0; i < nx; i++ ) 00574 { 00575 if ( ( ipts[i] = (PLINT *) malloc( (size_t) ny * sizeof ( PLINT * ) ) ) == NULL ) 00576 { 00577 plexit( "plfcont: Insufficient memory" ); 00578 } 00579 } 00580 00581 for ( i = 0; i < nlevel; i++ ) 00582 { 00583 plcntr( f2eval, f2eval_data, 00584 nx, ny, kx - 1, lx - 1, ky - 1, ly - 1, clevel[i], ipts, 00585 pltr, pltr_data ); 00586 00587 if ( error ) 00588 { 00589 error = 0; 00590 goto done; 00591 } 00592 } 00593 00594 done: 00595 for ( i = 0; i < nx; i++ ) 00596 { 00597 free( (void *) ipts[i] ); 00598 } 00599 free( (void *) ipts ); 00600 } 00601 00602 //-------------------------------------------------------------------------- 00603 // void plcntr() 00604 // 00605 // The contour for a given level is drawn here. Note iscan has nx 00606 // elements. ixstor and iystor each have nstor elements. 00607 //-------------------------------------------------------------------------- 00608 00609 static void 00610 plcntr( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ), 00611 PLPointer f2eval_data, 00612 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00613 PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts, 00614 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00615 PLPointer pltr_data ) 00616 { 00617 PLINT kcol, krow, lastindex; 00618 PLFLT distance; 00619 PLFLT save_def, save_scale; 00620 00621 char flabel[30]; 00622 plgchr( &save_def, &save_scale ); 00623 save_scale = save_scale / save_def; 00624 00625 cont_new_store( flev ); 00626 00627 // format contour label for plptex and define the font height of the labels 00628 plfloatlabel( flev, flabel, 30 ); 00629 plschr( 0.0, contlabel_size ); 00630 00631 // Clear array for traversed squares 00632 for ( kcol = kx; kcol < lx; kcol++ ) 00633 { 00634 for ( krow = ky; krow < ly; krow++ ) 00635 { 00636 ipts[kcol][krow] = 0; 00637 } 00638 } 00639 00640 00641 for ( krow = ky; krow < ly; krow++ ) 00642 { 00643 for ( kcol = kx; kcol < lx; kcol++ ) 00644 { 00645 if ( ipts[kcol][krow] == 0 ) 00646 { 00647 // Follow and draw a contour 00648 pldrawcn( f2eval, f2eval_data, 00649 nx, ny, kx, lx, ky, ly, flev, flabel, kcol, krow, 00650 0.0, 0.0, -2, ipts, &distance, &lastindex, 00651 pltr, pltr_data ); 00652 00653 if ( error ) 00654 return; 00655 } 00656 } 00657 } 00658 plschr( save_def, save_scale ); 00659 } 00660 00661 //-------------------------------------------------------------------------- 00662 // void pldrawcn() 00663 // 00664 // Follow and draw a contour. 00665 //-------------------------------------------------------------------------- 00666 00667 static void 00668 pldrawcn( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ), 00669 PLPointer f2eval_data, 00670 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00671 PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow, 00672 PLFLT lastx, PLFLT lasty, PLINT startedge, PLINT **ipts, 00673 PLFLT *distance, PLINT *lastindex, 00674 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00675 PLPointer pltr_data ) 00676 { 00677 PLFLT f[4]; 00678 PLFLT px[4], py[4], locx[4], locy[4]; 00679 PLINT iedge[4]; 00680 PLINT i, j, k, num, first, inext, kcolnext, krownext; 00681 00682 00683 ( *pltr )( kcol, krow + 1, &px[0], &py[0], pltr_data ); 00684 ( *pltr )( kcol, krow, &px[1], &py[1], pltr_data ); 00685 ( *pltr )( kcol + 1, krow, &px[2], &py[2], pltr_data ); 00686 ( *pltr )( kcol + 1, krow + 1, &px[3], &py[3], pltr_data ); 00687 00688 f[0] = f2eval( kcol, krow + 1, f2eval_data ) - flev; 00689 f[1] = f2eval( kcol, krow, f2eval_data ) - flev; 00690 f[2] = f2eval( kcol + 1, krow, f2eval_data ) - flev; 00691 f[3] = f2eval( kcol + 1, krow + 1, f2eval_data ) - flev; 00692 00693 for ( i = 0, j = 1; i < 4; i++, j = ( j + 1 ) % 4 ) 00694 { 00695 iedge[i] = ( f[i] * f[j] > 0.0 ) ? -1 : ( ( f[i] * f[j] < 0.0 ) ? 1 : 0 ); 00696 } 00697 00698 // Mark this square as done 00699 ipts[kcol][krow] = 1; 00700 00701 // Check if no contour has been crossed i.e. iedge[i] = -1 00702 if ( ( iedge[0] == -1 ) && ( iedge[1] == -1 ) && ( iedge[2] == -1 ) 00703 && ( iedge[3] == -1 ) ) 00704 return; 00705 00706 // Check if this is a completely flat square - in which case 00707 // ignore it 00708 if ( ( f[0] == 0.0 ) && ( f[1] == 0.0 ) && ( f[2] == 0.0 ) && 00709 ( f[3] == 0.0 ) ) 00710 return; 00711 00712 // Calculate intersection points 00713 num = 0; 00714 if ( startedge < 0 ) 00715 { 00716 first = 1; 00717 } 00718 else 00719 { 00720 locx[num] = lastx; 00721 locy[num] = lasty; 00722 num++; 00723 first = 0; 00724 } 00725 for ( k = 0, i = ( startedge < 0 ? 0 : startedge ); k < 4; k++, i = ( i + 1 ) % 4 ) 00726 { 00727 if ( i == startedge ) 00728 continue; 00729 00730 // If the contour is an edge check it hasn't already been done 00731 if ( f[i] == 0.0 && f[( i + 1 ) % 4] == 0.0 ) 00732 { 00733 kcolnext = kcol; 00734 krownext = krow; 00735 if ( i == 0 ) 00736 kcolnext--; 00737 if ( i == 1 ) 00738 krownext--; 00739 if ( i == 2 ) 00740 kcolnext++; 00741 if ( i == 3 ) 00742 krownext++; 00743 if ( ( kcolnext < kx ) || ( kcolnext >= lx ) || 00744 ( krownext < ky ) || ( krownext >= ly ) || 00745 ( ipts[kcolnext][krownext] == 1 ) ) 00746 continue; 00747 } 00748 if ( ( iedge[i] == 1 ) || ( f[i] == 0.0 ) ) 00749 { 00750 j = ( i + 1 ) % 4; 00751 if ( f[i] != 0.0 ) 00752 { 00753 locx[num] = ( px[i] * fabs( f[j] ) + px[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] ); 00754 locy[num] = ( py[i] * fabs( f[j] ) + py[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] ); 00755 } 00756 else 00757 { 00758 locx[num] = px[i]; 00759 locy[num] = py[i]; 00760 } 00761 // If this is the start of the contour then move to the point 00762 if ( first == 1 ) 00763 { 00764 cont_mv_store( locx[num], locy[num] ); 00765 first = 0; 00766 *distance = 0; 00767 *lastindex = 0; 00768 } 00769 else 00770 { 00771 // Link to the next point on the contour 00772 if ( contlabel_active ) 00773 pl_drawcontlabel( locx[num], locy[num], flabel, distance, lastindex ); 00774 else 00775 cont_xy_store( locx[num], locy[num] ); 00776 // Need to follow contour into next grid box 00777 // Easy case where contour does not pass through corner 00778 if ( f[i] != 0.0 ) 00779 { 00780 kcolnext = kcol; 00781 krownext = krow; 00782 inext = ( i + 2 ) % 4; 00783 if ( i == 0 ) 00784 kcolnext--; 00785 if ( i == 1 ) 00786 krownext--; 00787 if ( i == 2 ) 00788 kcolnext++; 00789 if ( i == 3 ) 00790 krownext++; 00791 if ( ( kcolnext >= kx ) && ( kcolnext < lx ) && 00792 ( krownext >= ky ) && ( krownext < ly ) && 00793 ( ipts[kcolnext][krownext] == 0 ) ) 00794 { 00795 pldrawcn( f2eval, f2eval_data, 00796 nx, ny, kx, lx, ky, ly, flev, flabel, 00797 kcolnext, krownext, 00798 locx[num], locy[num], inext, ipts, 00799 distance, lastindex, 00800 pltr, pltr_data ); 00801 } 00802 } 00803 // Hard case where contour passes through corner 00804 // This is still not perfect - it may lose the contour 00805 // which won't upset the contour itself (we can find it 00806 // again later) but might upset the labelling 00807 else 00808 { 00809 kcolnext = kcol; 00810 krownext = krow; 00811 inext = ( i + 2 ) % 4; 00812 if ( i == 0 ) 00813 { 00814 kcolnext--; krownext++; 00815 } 00816 if ( i == 1 ) 00817 { 00818 krownext--; kcolnext--; 00819 } 00820 if ( i == 2 ) 00821 { 00822 kcolnext++; krownext--; 00823 } 00824 if ( i == 3 ) 00825 { 00826 krownext++; kcolnext++; 00827 } 00828 if ( ( kcolnext >= kx ) && ( kcolnext < lx ) && 00829 ( krownext >= ky ) && ( krownext < ly ) && 00830 ( ipts[kcolnext][krownext] == 0 ) ) 00831 { 00832 pldrawcn( f2eval, f2eval_data, 00833 nx, ny, kx, lx, ky, ly, flev, flabel, 00834 kcolnext, krownext, 00835 locx[num], locy[num], inext, ipts, 00836 distance, lastindex, 00837 pltr, pltr_data ); 00838 } 00839 } 00840 if ( first == 1 ) 00841 { 00842 // Move back to first point 00843 cont_mv_store( locx[num], locy[num] ); 00844 first = 0; 00845 *distance = 0; 00846 *lastindex = 0; 00847 first = 0; 00848 } 00849 else 00850 { 00851 first = 1; 00852 } 00853 num++; 00854 } 00855 } 00856 } 00857 } 00858 00859 //-------------------------------------------------------------------------- 00860 // pltr0() 00861 // 00862 // Identity transformation. 00863 //-------------------------------------------------------------------------- 00864 00865 void 00866 pltr0( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer PL_UNUSED( pltr_data ) ) 00867 { 00868 *tx = x; 00869 *ty = y; 00870 } 00871 00872 //-------------------------------------------------------------------------- 00873 // pltr1() 00874 // 00875 // Does linear interpolation from singly dimensioned coord arrays. 00876 // 00877 // Just abort for now if coordinates are out of bounds (don't think it's 00878 // possible, but if so we could use linear extrapolation). 00879 //-------------------------------------------------------------------------- 00880 00881 void 00882 pltr1( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) 00883 { 00884 PLINT ul, ur, vl, vr; 00885 PLFLT du, dv; 00886 PLFLT xl, xr, yl, yr; 00887 00888 PLcGrid *grid = (PLcGrid *) pltr_data; 00889 PLFLT *xg = grid->xg; 00890 PLFLT *yg = grid->yg; 00891 PLINT nx = grid->nx; 00892 PLINT ny = grid->ny; 00893 00894 ul = (PLINT) x; 00895 ur = ul + 1; 00896 du = x - ul; 00897 00898 vl = (PLINT) y; 00899 vr = vl + 1; 00900 dv = y - vl; 00901 00902 if ( x < 0 || x > nx - 1 || y < 0 || y > ny - 1 ) 00903 { 00904 plexit( "pltr1: Invalid coordinates" ); 00905 } 00906 00907 // Look up coordinates in row-dominant array. 00908 // Have to handle right boundary specially -- if at the edge, we'd better 00909 // not reference the out of bounds point. 00910 // 00911 00912 xl = xg[ul]; 00913 yl = yg[vl]; 00914 00915 if ( ur == nx ) 00916 { 00917 *tx = xl; 00918 } 00919 else 00920 { 00921 xr = xg[ur]; 00922 *tx = xl * ( 1 - du ) + xr * du; 00923 } 00924 if ( vr == ny ) 00925 { 00926 *ty = yl; 00927 } 00928 else 00929 { 00930 yr = yg[vr]; 00931 *ty = yl * ( 1 - dv ) + yr * dv; 00932 } 00933 } 00934 00935 //-------------------------------------------------------------------------- 00936 // pltr2() 00937 // 00938 // Does linear interpolation from doubly dimensioned coord arrays (column 00939 // dominant, as per normal C 2d arrays). 00940 // 00941 // This routine includes lots of checks for out of bounds. This would occur 00942 // occasionally due to some bugs in the contour plotter (now fixed). If an 00943 // out of bounds coordinate is obtained, the boundary value is provided 00944 // along with a warning. These checks should stay since no harm is done if 00945 // if everything works correctly. 00946 //-------------------------------------------------------------------------- 00947 00948 void 00949 pltr2( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) 00950 { 00951 PLINT ul, ur, vl, vr; 00952 PLFLT du, dv; 00953 PLFLT xll, xlr, xrl, xrr; 00954 PLFLT yll, ylr, yrl, yrr; 00955 PLFLT xmin, xmax, ymin, ymax; 00956 00957 PLcGrid2 *grid = (PLcGrid2 *) pltr_data; 00958 PLFLT **xg = grid->xg; 00959 PLFLT **yg = grid->yg; 00960 PLINT nx = grid->nx; 00961 PLINT ny = grid->ny; 00962 00963 ul = (PLINT) x; 00964 ur = ul + 1; 00965 du = x - ul; 00966 00967 vl = (PLINT) y; 00968 vr = vl + 1; 00969 dv = y - vl; 00970 00971 xmin = 0; 00972 xmax = nx - 1; 00973 ymin = 0; 00974 ymax = ny - 1; 00975 00976 if ( x < xmin || x > xmax || y < ymin || y > ymax ) 00977 { 00978 plwarn( "pltr2: Invalid coordinates" ); 00979 if ( x < xmin ) 00980 { 00981 if ( y < ymin ) 00982 { 00983 *tx = xg[0][0]; 00984 *ty = yg[0][0]; 00985 } 00986 else if ( y > ymax ) 00987 { 00988 *tx = xg[0][ny - 1]; 00989 *ty = yg[0][ny - 1]; 00990 } 00991 else 00992 { 00993 xll = xg[0][vl]; 00994 yll = yg[0][vl]; 00995 xlr = xg[0][vr]; 00996 ylr = yg[0][vr]; 00997 00998 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 00999 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01000 } 01001 } 01002 else if ( x > xmax ) 01003 { 01004 if ( y < ymin ) 01005 { 01006 *tx = xg[nx - 1][0]; 01007 *ty = yg[nx - 1][0]; 01008 } 01009 else if ( y > ymax ) 01010 { 01011 *tx = xg[nx - 1][ny - 1]; 01012 *ty = yg[nx - 1][ny - 1]; 01013 } 01014 else 01015 { 01016 xll = xg[nx - 1][vl]; 01017 yll = yg[nx - 1][vl]; 01018 xlr = xg[nx - 1][vr]; 01019 ylr = yg[nx - 1][vr]; 01020 01021 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01022 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01023 } 01024 } 01025 else 01026 { 01027 if ( y < ymin ) 01028 { 01029 xll = xg[ul][0]; 01030 xrl = xg[ur][0]; 01031 yll = yg[ul][0]; 01032 yrl = yg[ur][0]; 01033 01034 *tx = xll * ( 1 - du ) + xrl * ( du ); 01035 *ty = yll * ( 1 - du ) + yrl * ( du ); 01036 } 01037 else if ( y > ymax ) 01038 { 01039 xlr = xg[ul][ny - 1]; 01040 xrr = xg[ur][ny - 1]; 01041 ylr = yg[ul][ny - 1]; 01042 yrr = yg[ur][ny - 1]; 01043 01044 *tx = xlr * ( 1 - du ) + xrr * ( du ); 01045 *ty = ylr * ( 1 - du ) + yrr * ( du ); 01046 } 01047 } 01048 } 01049 01050 // Normal case. 01051 // Look up coordinates in row-dominant array. 01052 // Have to handle right boundary specially -- if at the edge, we'd 01053 // better not reference the out of bounds point. 01054 // 01055 01056 else 01057 { 01058 xll = xg[ul][vl]; 01059 yll = yg[ul][vl]; 01060 01061 // ur is out of bounds 01062 01063 if ( ur == nx && vr < ny ) 01064 { 01065 xlr = xg[ul][vr]; 01066 ylr = yg[ul][vr]; 01067 01068 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01069 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01070 } 01071 01072 // vr is out of bounds 01073 01074 else if ( ur < nx && vr == ny ) 01075 { 01076 xrl = xg[ur][vl]; 01077 yrl = yg[ur][vl]; 01078 01079 *tx = xll * ( 1 - du ) + xrl * ( du ); 01080 *ty = yll * ( 1 - du ) + yrl * ( du ); 01081 } 01082 01083 // both ur and vr are out of bounds 01084 01085 else if ( ur == nx && vr == ny ) 01086 { 01087 *tx = xll; 01088 *ty = yll; 01089 } 01090 01091 // everything in bounds 01092 01093 else 01094 { 01095 xrl = xg[ur][vl]; 01096 xlr = xg[ul][vr]; 01097 xrr = xg[ur][vr]; 01098 01099 yrl = yg[ur][vl]; 01100 ylr = yg[ul][vr]; 01101 yrr = yg[ur][vr]; 01102 01103 *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + 01104 xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); 01105 01106 *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + 01107 yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); 01108 } 01109 } 01110 } 01111 01112 //-------------------------------------------------------------------------- 01113 // pltr2p() 01114 // 01115 // Just like pltr2() but uses pointer arithmetic to get coordinates from 2d 01116 // grid tables. This form of grid tables is compatible with those from 01117 // PLplot 4.0. The grid data must be pointed to by a PLcGrid structure. 01118 //-------------------------------------------------------------------------- 01119 01120 void 01121 pltr2p( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) 01122 { 01123 PLINT ul, ur, vl, vr; 01124 PLFLT du, dv; 01125 PLFLT xll, xlr, xrl, xrr; 01126 PLFLT yll, ylr, yrl, yrr; 01127 PLFLT xmin, xmax, ymin, ymax; 01128 01129 PLcGrid *grid = (PLcGrid *) pltr_data; 01130 PLFLT *xg = grid->xg; 01131 PLFLT *yg = grid->yg; 01132 PLINT nx = grid->nx; 01133 PLINT ny = grid->ny; 01134 01135 ul = (PLINT) x; 01136 ur = ul + 1; 01137 du = x - ul; 01138 01139 vl = (PLINT) y; 01140 vr = vl + 1; 01141 dv = y - vl; 01142 01143 xmin = 0; 01144 xmax = nx - 1; 01145 ymin = 0; 01146 ymax = ny - 1; 01147 01148 if ( x < xmin || x > xmax || y < ymin || y > ymax ) 01149 { 01150 plwarn( "pltr2p: Invalid coordinates" ); 01151 if ( x < xmin ) 01152 { 01153 if ( y < ymin ) 01154 { 01155 *tx = *xg; 01156 *ty = *yg; 01157 } 01158 else if ( y > ymax ) 01159 { 01160 *tx = *( xg + ( ny - 1 ) ); 01161 *ty = *( yg + ( ny - 1 ) ); 01162 } 01163 else 01164 { 01165 ul = 0; 01166 xll = *( xg + ul * ny + vl ); 01167 yll = *( yg + ul * ny + vl ); 01168 xlr = *( xg + ul * ny + vr ); 01169 ylr = *( yg + ul * ny + vr ); 01170 01171 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01172 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01173 } 01174 } 01175 else if ( x > xmax ) 01176 { 01177 if ( y < ymin ) 01178 { 01179 *tx = *( xg + ( ny - 1 ) * nx ); 01180 *ty = *( yg + ( ny - 1 ) * nx ); 01181 } 01182 else if ( y > ymax ) 01183 { 01184 *tx = *( xg + ( ny - 1 ) + ( nx - 1 ) * ny ); 01185 *ty = *( yg + ( ny - 1 ) + ( nx - 1 ) * ny ); 01186 } 01187 else 01188 { 01189 ul = nx - 1; 01190 xll = *( xg + ul * ny + vl ); 01191 yll = *( yg + ul * ny + vl ); 01192 xlr = *( xg + ul * ny + vr ); 01193 ylr = *( yg + ul * ny + vr ); 01194 01195 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01196 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01197 } 01198 } 01199 else 01200 { 01201 if ( y < ymin ) 01202 { 01203 vl = 0; 01204 xll = *( xg + ul * ny + vl ); 01205 xrl = *( xg + ur * ny + vl ); 01206 yll = *( yg + ul * ny + vl ); 01207 yrl = *( yg + ur * ny + vl ); 01208 01209 *tx = xll * ( 1 - du ) + xrl * ( du ); 01210 *ty = yll * ( 1 - du ) + yrl * ( du ); 01211 } 01212 else if ( y > ymax ) 01213 { 01214 vr = ny - 1; 01215 xlr = *( xg + ul * ny + vr ); 01216 xrr = *( xg + ur * ny + vr ); 01217 ylr = *( yg + ul * ny + vr ); 01218 yrr = *( yg + ur * ny + vr ); 01219 01220 *tx = xlr * ( 1 - du ) + xrr * ( du ); 01221 *ty = ylr * ( 1 - du ) + yrr * ( du ); 01222 } 01223 } 01224 } 01225 01226 // Normal case. 01227 // Look up coordinates in row-dominant array. 01228 // Have to handle right boundary specially -- if at the edge, we'd better 01229 // not reference the out of bounds point. 01230 // 01231 01232 else 01233 { 01234 xll = *( xg + ul * ny + vl ); 01235 yll = *( yg + ul * ny + vl ); 01236 01237 // ur is out of bounds 01238 01239 if ( ur == nx && vr < ny ) 01240 { 01241 xlr = *( xg + ul * ny + vr ); 01242 ylr = *( yg + ul * ny + vr ); 01243 01244 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01245 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01246 } 01247 01248 // vr is out of bounds 01249 01250 else if ( ur < nx && vr == ny ) 01251 { 01252 xrl = *( xg + ur * ny + vl ); 01253 yrl = *( yg + ur * ny + vl ); 01254 01255 *tx = xll * ( 1 - du ) + xrl * ( du ); 01256 *ty = yll * ( 1 - du ) + yrl * ( du ); 01257 } 01258 01259 // both ur and vr are out of bounds 01260 01261 else if ( ur == nx && vr == ny ) 01262 { 01263 *tx = xll; 01264 *ty = yll; 01265 } 01266 01267 // everything in bounds 01268 01269 else 01270 { 01271 xrl = *( xg + ur * ny + vl ); 01272 xlr = *( xg + ul * ny + vr ); 01273 xrr = *( xg + ur * ny + vr ); 01274 01275 yrl = *( yg + ur * ny + vl ); 01276 ylr = *( yg + ul * ny + vr ); 01277 yrr = *( yg + ur * ny + vr ); 01278 01279 *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + 01280 xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); 01281 01282 *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + 01283 yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); 01284 } 01285 } 01286 } 01287 01288 //-------------------------------------------------------------------------- 01289 // pltr2f() 01290 // 01291 // Does linear interpolation from doubly dimensioned coord arrays 01292 // (row dominant, i.e. Fortran ordering). 01293 // 01294 // This routine includes lots of checks for out of bounds. This would 01295 // occur occasionally due to a bug in the contour plotter that is now fixed. 01296 // If an out of bounds coordinate is obtained, the boundary value is provided 01297 // along with a warning. These checks should stay since no harm is done if 01298 // if everything works correctly. 01299 //-------------------------------------------------------------------------- 01300 01301 void 01302 pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data ) 01303 { 01304 PLINT ul, ur, vl, vr; 01305 PLFLT du, dv; 01306 PLFLT xll, xlr, xrl, xrr; 01307 PLFLT yll, ylr, yrl, yrr; 01308 PLFLT xmin, xmax, ymin, ymax; 01309 01310 PLcGrid *cgrid = (PLcGrid *) pltr_data; 01311 PLFLT *xg = cgrid->xg; 01312 PLFLT *yg = cgrid->yg; 01313 PLINT nx = cgrid->nx; 01314 PLINT ny = cgrid->ny; 01315 01316 ul = (PLINT) x; 01317 ur = ul + 1; 01318 du = x - ul; 01319 01320 vl = (PLINT) y; 01321 vr = vl + 1; 01322 dv = y - vl; 01323 01324 xmin = 0; 01325 xmax = nx - 1; 01326 ymin = 0; 01327 ymax = ny - 1; 01328 01329 if ( x < xmin || x > xmax || y < ymin || y > ymax ) 01330 { 01331 plwarn( "pltr2f: Invalid coordinates" ); 01332 01333 if ( x < xmin ) 01334 { 01335 if ( y < ymin ) 01336 { 01337 *tx = *xg; 01338 *ty = *yg; 01339 } 01340 else if ( y > ymax ) 01341 { 01342 *tx = *( xg + ( ny - 1 ) * nx ); 01343 *ty = *( yg + ( ny - 1 ) * nx ); 01344 } 01345 else 01346 { 01347 ul = 0; 01348 xll = *( xg + ul + vl * nx ); 01349 yll = *( yg + ul + vl * nx ); 01350 xlr = *( xg + ul + vr * nx ); 01351 ylr = *( yg + ul + vr * nx ); 01352 01353 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01354 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01355 } 01356 } 01357 else if ( x > xmax ) 01358 { 01359 if ( y < ymin ) 01360 { 01361 *tx = *( xg + ( nx - 1 ) ); 01362 *ty = *( yg + ( nx - 1 ) ); 01363 } 01364 else if ( y > ymax ) 01365 { 01366 *tx = *( xg + ( nx - 1 ) + ( ny - 1 ) * nx ); 01367 *ty = *( yg + ( nx - 1 ) + ( ny - 1 ) * nx ); 01368 } 01369 else 01370 { 01371 ul = nx - 1; 01372 xll = *( xg + ul + vl * nx ); 01373 yll = *( yg + ul + vl * nx ); 01374 xlr = *( xg + ul + vr * nx ); 01375 ylr = *( yg + ul + vr * nx ); 01376 01377 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01378 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01379 } 01380 } 01381 else 01382 { 01383 if ( y < ymin ) 01384 { 01385 vl = 0; 01386 xll = *( xg + ul + vl * nx ); 01387 xrl = *( xg + ur + vl * nx ); 01388 yll = *( yg + ul + vl * nx ); 01389 yrl = *( yg + ur + vl * nx ); 01390 01391 *tx = xll * ( 1 - du ) + xrl * ( du ); 01392 *ty = yll * ( 1 - du ) + yrl * ( du ); 01393 } 01394 else if ( y > ymax ) 01395 { 01396 vr = ny - 1; 01397 xlr = *( xg + ul + vr * nx ); 01398 xrr = *( xg + ur + vr * nx ); 01399 ylr = *( yg + ul + vr * nx ); 01400 yrr = *( yg + ur + vr * nx ); 01401 01402 *tx = xlr * ( 1 - du ) + xrr * ( du ); 01403 *ty = ylr * ( 1 - du ) + yrr * ( du ); 01404 } 01405 } 01406 } 01407 01408 // Normal case. 01409 // Look up coordinates in row-dominant array. 01410 // Have to handle right boundary specially -- if at the edge, we'd 01411 // better not reference the out of bounds point. 01412 01413 else 01414 { 01415 xll = *( xg + ul + vl * nx ); 01416 yll = *( yg + ul + vl * nx ); 01417 01418 // ur is out of bounds 01419 01420 if ( ur == nx && vr < ny ) 01421 { 01422 xlr = *( xg + ul + vr * nx ); 01423 ylr = *( yg + ul + vr * nx ); 01424 01425 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01426 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01427 } 01428 01429 // vr is out of bounds 01430 01431 else if ( ur < nx && vr == ny ) 01432 { 01433 xrl = *( xg + ur + vl * nx ); 01434 yrl = *( yg + ur + vl * nx ); 01435 01436 *tx = xll * ( 1 - du ) + xrl * ( du ); 01437 *ty = yll * ( 1 - du ) + yrl * ( du ); 01438 } 01439 01440 // both ur and vr are out of bounds 01441 01442 else if ( ur == nx && vr == ny ) 01443 { 01444 *tx = xll; 01445 *ty = yll; 01446 } 01447 01448 // everything in bounds 01449 01450 else 01451 { 01452 xrl = *( xg + ur + vl * nx ); 01453 xlr = *( xg + ul + vr * nx ); 01454 xrr = *( xg + ur + vr * nx ); 01455 01456 yrl = *( yg + ur + vl * nx ); 01457 ylr = *( yg + ul + vr * nx ); 01458 yrr = *( yg + ur + vr * nx ); 01459 *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + 01460 xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); 01461 01462 *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + 01463 yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); 01464 } 01465 } 01466 }