PLplot  5.10.0
plcont.c
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines