PLplot  5.11.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
plcont.c
Go to the documentation of this file.
1 // Contour plotter.
2 //
3 // Copyright (C) 1995, 2000, 2001 Maurice LeBrun
4 // Copyright (C) 2000, 2002 Joao Cardoso
5 // Copyright (C) 2000-2014 Alan W. Irwin
6 // Copyright (C) 2004 Andrew Ross
7 //
8 // This file is part of PLplot.
9 //
10 // PLplot is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Library General Public License as published
12 // by the Free Software Foundation; either version 2 of the License, or
13 // (at your option) any later version.
14 //
15 // PLplot is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Library General Public License for more details.
19 //
20 // You should have received a copy of the GNU Library General Public License
21 // along with PLplot; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 //
24 
25 #include "plplotP.h"
26 
27 #ifdef MSDOS
28 #pragma optimize("",off)
29 #endif
30 
31 // Static function prototypes.
32 
33 static void
35  PLPointer plf2eval_data,
36  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
37  PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts,
38  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
39  PLPointer pltr_data );
40 
41 static void
43  PLPointer plf2eval_data,
44  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
45  PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow,
46  PLFLT lastx, PLFLT lasty, PLINT startedge,
47  PLINT **ipts, PLFLT *distance, PLINT *lastindex,
48  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
49  PLPointer pltr_data );
50 
51 static void
52 plfloatlabel( PLFLT value, char *string, PLINT len );
53 
54 static PLFLT
55 plP_pcwcx( PLINT x );
56 
57 static PLFLT
58 plP_pcwcy( PLINT y );
59 
60 static void
61 pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex );
62 
63 // Error flag for aborts
64 
65 static int error;
66 
67 //**************************************
68 //
69 // Defaults for contour label printing.
70 //
71 //**************************************
72 
73 // Font height for contour labels (normalized)
74 static PLFLT
76 
77 // Offset of label from contour line (if set to 0.0, labels are printed on the lines).
78 static PLFLT
80 
81 // Spacing parameter for contour labels
82 static PLFLT
84 
85 // Activate labels, default off
86 static PLINT
88 
89 // If the contour label exceed 10^(limexp) or 10^(-limexp), the exponential format is used
90 static PLINT
91  limexp = 4;
92 
93 // Number of significant digits
94 static PLINT
95  sigprec = 2;
96 
97 //******* contour lines storage ***************************
98 
99 static CONT_LEVEL *startlev = NULL;
102 
103 static int cont3d = 0;
104 
105 static CONT_LINE *
106 alloc_line( void )
107 {
108  CONT_LINE *line;
109 
110  if ( ( line = (CONT_LINE *) malloc( sizeof ( CONT_LINE ) ) ) == NULL )
111  {
112  plexit( "alloc_line: Insufficient memory" );
113  }
114 
115  line->x = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) );
116  line->y = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) );
117 
118  if ( ( line->x == NULL ) || ( line->y == NULL ) )
119  {
120  plexit( "alloc_line: Insufficient memory" );
121  }
122 
123  line->npts = 0;
124  line->next = NULL;
125 
126  return line;
127 }
128 
129 static CONT_LEVEL *
131 {
132  CONT_LEVEL *node;
133 
134  if ( ( node = (CONT_LEVEL *) malloc( sizeof ( CONT_LEVEL ) ) ) == NULL )
135  {
136  plexit( "alloc_level: Insufficient memory" );
137  }
138  node->level = level;
139  node->next = NULL;
140  node->line = alloc_line( );
141 
142  return node;
143 }
144 
145 static void
147 {
148  if ( ( ( line->x = (PLFLT *) realloc( line->x,
149  (size_t) ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) ||
150  ( ( line->y = (PLFLT *) realloc( line->y,
151  (size_t) ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) )
152  plexit( "realloc_line: Insufficient memory" );
153 }
154 
155 
156 // new contour level
157 static void
159 {
160  if ( cont3d )
161  {
162  if ( startlev == NULL )
163  {
164  startlev = alloc_level( level );
165  currlev = startlev;
166  }
167  else
168  {
169  currlev->next = alloc_level( level );
170  currlev = currlev->next;
171  }
172  currline = currlev->line;
173  }
174 }
175 
176 void
178 {
179  CONT_LINE *tline, *cline;
180  CONT_LEVEL *tlev, *clevel;
181 
182  if ( ct != NULL )
183  {
184  clevel = ct;
185 
186  do
187  {
188  cline = clevel->line;
189  do
190  {
191 #ifdef CONT_PLOT_DEBUG // for 2D plots. For 3D plots look at plot3.c:plotsh3di()
192  plP_movwor( cline->x[0], cline->y[0] );
193  for ( j = 1; j < cline->npts; j++ )
194  plP_drawor( cline->x[j], cline->y[j] );
195 #endif
196  tline = cline->next;
197  free( cline->x );
198  free( cline->y );
199  free( cline );
200  cline = tline;
201  }
202  while ( cline != NULL );
203  tlev = clevel->next;
204  free( clevel );
205  clevel = tlev;
206  }
207  while ( clevel != NULL );
208  startlev = NULL;
209  }
210 }
211 
212 static void
214 {
215  if ( cont3d )
216  {
217  PLINT pts = currline->npts;
218 
219  if ( pts % LINE_ITEMS == 0 )
220  realloc_line( currline );
221 
222  currline->x[pts] = xx;
223  currline->y[pts] = yy;
224  currline->npts++;
225  }
226  else
227  plP_drawor( xx, yy );
228 }
229 
230 static void
232 {
233  if ( cont3d )
234  {
235  if ( currline->npts != 0 ) // not an empty list, allocate new
236  {
237  currline->next = alloc_line( );
238  currline = currline->next;
239  }
240 
241  // and fill first element
242  currline->x[0] = xx;
243  currline->y[0] = yy;
244  currline->npts = 1;
245  }
246  else
247  plP_movwor( xx, yy );
248 }
249 
250 // small routine to set offset and spacing of contour labels, see desciption above
251 void c_pl_setcontlabelparam( PLFLT offset, PLFLT size, PLFLT spacing, PLINT active )
252 {
253  contlabel_offset = offset;
254  contlabel_size = size;
255  contlabel_space = spacing;
256  contlabel_active = active;
257 }
258 
259 // small routine to set the format of the contour labels, description of limexp and prec see above
260 void c_pl_setcontlabelformat( PLINT lexp, PLINT sigdig )
261 {
262  limexp = lexp;
263  sigprec = sigdig;
264 }
265 
266 static void pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex )
267 {
268  PLFLT delta_x, delta_y;
269  PLINT currx_old, curry_old;
270 
271  delta_x = plP_pcdcx( plsc->currx ) - plP_pcdcx( plP_wcpcx( tpx ) );
272  delta_y = plP_pcdcy( plsc->curry ) - plP_pcdcy( plP_wcpcy( tpy ) );
273 
274  currx_old = plsc->currx;
275  curry_old = plsc->curry;
276 
277  *distance += sqrt( delta_x * delta_x + delta_y * delta_y );
278 
279  plP_drawor( tpx, tpy );
280 
281  if ( (int) ( fabs( *distance / contlabel_space ) ) > *lastindex )
282  {
283  PLFLT scale, vec_x, vec_y, mx, my, dev_x, dev_y, off_x, off_y;
284 
285  vec_x = tpx - plP_pcwcx( currx_old );
286  vec_y = tpy - plP_pcwcy( curry_old );
287 
288  // Ensure labels appear the right way up
289  if ( vec_x < 0 )
290  {
291  vec_x = -vec_x;
292  vec_y = -vec_y;
293  }
294 
295  mx = (double) plsc->wpxscl / (double) plsc->phyxlen;
296  my = (double) plsc->wpyscl / (double) plsc->phyylen;
297 
298  dev_x = -my * vec_y / mx;
299  dev_y = mx * vec_x / my;
300 
301  scale = sqrt( ( mx * mx * dev_x * dev_x + my * my * dev_y * dev_y ) /
303 
304  off_x = dev_x / scale;
305  off_y = dev_y / scale;
306 
307  plptex( tpx + off_x, tpy + off_y, vec_x, vec_y, 0.5, flabel );
308  plP_movwor( tpx, tpy );
309  ( *lastindex )++;
310  }
311  else
312  plP_movwor( tpx, tpy );
313 }
314 
315 
316 // Format contour labels. Arguments:
317 // value: floating point number to be formatted
318 // string: the formatted label, plptex must be called with it to actually
319 // print the label
320 //
321 
322 static void plfloatlabel( PLFLT value, char *string, PLINT len )
323 {
324  PLINT setpre, precis;
325  // form[10] gives enough space for all non-malicious formats.
326  // tmpstring[15] gives enough room for 3 digits in a negative exponent
327  // or 4 digits in a positive exponent + null termination. That
328  // should be enough for all non-malicious use.
329  // Obviously there are security issues here that
330  // should be addressed as well.
331  //
332 #define FORM_LEN 10
333 #define TMPSTRING_LEN 15
334  char form[FORM_LEN], tmpstring[TMPSTRING_LEN];
335  PLINT exponent = 0;
336  PLFLT mant, tmp;
337 
338  PLINT prec = sigprec;
339 
340  plP_gprec( &setpre, &precis );
341 
342  if ( setpre )
343  prec = precis;
344 
345  if ( value > 0.0 )
346  tmp = log10( value );
347  else if ( value < 0.0 )
348  tmp = log10( -value );
349  else
350  tmp = 0;
351 
352  if ( tmp >= 0.0 )
353  exponent = (int) tmp;
354  else if ( tmp < 0.0 )
355  {
356  tmp = -tmp;
357  if ( floor( tmp ) < tmp )
358  exponent = -(int) ( floor( tmp ) + 1.0 );
359  else
360  exponent = -(int) ( floor( tmp ) );
361  }
362 
363  mant = value / pow( 10.0, exponent );
364 
365  if ( mant != 0.0 )
366  mant = (int) ( mant * pow( 10.0, prec - 1 ) + 0.5 * mant / fabs( mant ) ) / pow( 10.0, prec - 1 );
367 
368  snprintf( form, FORM_LEN, "%%.%df", prec - 1 );
369  snprintf( string, (size_t) len, form, mant );
370  snprintf( tmpstring, TMPSTRING_LEN, "#(229)10#u%d", exponent );
371  strncat( string, tmpstring, (size_t) len - strlen( string ) - 1 );
372 
373  if ( abs( exponent ) < limexp || value == 0.0 )
374  {
375  value = pow( 10.0, exponent ) * mant;
376 
377  if ( exponent >= 0 )
378  prec = prec - 1 - exponent;
379  else
380  prec = prec - 1 + abs( exponent );
381 
382  if ( prec < 0 )
383  prec = 0;
384 
385  snprintf( form, FORM_LEN, "%%.%df", (int) prec );
386  snprintf( string, (size_t) len, form, value );
387  }
388 }
389 
390 // physical coords (x) to world coords
391 
392 static PLFLT
394 {
395  return ( ( x - plsc->wpxoff ) / plsc->wpxscl );
396 }
397 
398 // physical coords (y) to world coords
399 
400 static PLFLT
402 {
403  return ( ( y - plsc->wpyoff ) / plsc->wpyscl );
404 }
405 
406 //--------------------------------------------------------------------------
407 // plf2eval1()
408 //
409 // Does a lookup from a 2d function array. Array is of type (PLFLT **),
410 // and is column dominant (normal C ordering).
411 //--------------------------------------------------------------------------
412 
413 PLFLT
414 plf2eval1( PLINT ix, PLINT iy, PLPointer plf2eval_data )
415 {
416  PLFLT value;
417  PLFLT **z = (PLFLT **) plf2eval_data;
418 
419  value = z[ix][iy];
420 
421  return value;
422 }
423 
424 //--------------------------------------------------------------------------
425 // plf2eval2()
426 //
427 // Does a lookup from a 2d function array. plf2eval_data is treated as type
428 // (PLfGrid2 *).
429 //--------------------------------------------------------------------------
430 
431 PLFLT
432 plf2eval2( PLINT ix, PLINT iy, PLPointer plf2eval_data )
433 {
434  PLFLT value;
435  PLfGrid2 *grid = (PLfGrid2 *) plf2eval_data;
436 
437  value = grid->f[ix][iy];
438 
439  return value;
440 }
441 
442 //--------------------------------------------------------------------------
443 // plf2eval()
444 //
445 // Does a lookup from a 2d function array. Array is of type (PLFLT *), and
446 // is column dominant (normal C ordering). You MUST fill the ny maximum
447 // array index entry in the PLfGrid struct.
448 //--------------------------------------------------------------------------
449 
450 PLFLT
451 plf2eval( PLINT ix, PLINT iy, PLPointer plf2eval_data )
452 {
453  PLFLT value;
454  PLfGrid *grid = (PLfGrid *) plf2eval_data;
455 
456  value = grid->f[ix * grid->ny + iy];
457 
458  return value;
459 }
460 
461 //--------------------------------------------------------------------------
462 // plf2evalr()
463 //
464 // Does a lookup from a 2d function array. Array is of type (PLFLT *), and
465 // is row dominant (Fortran ordering). You MUST fill the nx maximum array
466 // index entry in the PLfGrid struct.
467 //--------------------------------------------------------------------------
468 
469 PLFLT
470 plf2evalr( PLINT ix, PLINT iy, PLPointer plf2eval_data )
471 {
472  PLFLT value;
473  PLfGrid *grid = (PLfGrid *) plf2eval_data;
474 
475  value = grid->f[ix + iy * grid->nx];
476 
477  return value;
478 }
479 
480 //--------------------------------------------------------------------------
481 //
482 // cont_store:
483 //
484 // Draw contour lines in memory.
485 // cont_clean_store() must be called after use to release allocated memory.
486 //
487 //--------------------------------------------------------------------------
488 
489 void
490 cont_store( const PLFLT * const *f, PLINT nx, PLINT ny, PLINT kx, PLINT lx,
491  PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel,
492  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
493  PLPointer pltr_data,
494  CONT_LEVEL **contour )
495 {
496  cont3d = 1;
497 
498  plcont( f, nx, ny, kx, lx, ky, ly, clevel, nlevel,
499  pltr, pltr_data );
500 
501  *contour = startlev;
502  cont3d = 0;
503 }
504 
505 //--------------------------------------------------------------------------
506 // void plcont()
507 //
508 // Draws a contour plot from data in f(nx,ny). Is just a front-end to
509 // plfcont, with a particular choice for f2eval and f2eval_data.
510 //--------------------------------------------------------------------------
511 
512 void
513 c_plcont( const PLFLT * const *f, PLINT nx, PLINT ny, PLINT kx, PLINT lx,
514  PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel,
515  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
516  PLPointer pltr_data )
517 {
518  if ( pltr == NULL )
519  {
520  // If pltr is undefined, abort with an error.
521  plabort( "plcont: The pltr callback must be defined" );
522  return;
523  }
524 
526  nx, ny, kx, lx, ky, ly, clevel, nlevel,
527  pltr, pltr_data );
528 }
529 
530 //--------------------------------------------------------------------------
531 // void plfcont()
532 //
533 // Draws a contour plot using the function evaluator f2eval and data stored
534 // by way of the f2eval_data pointer. This allows arbitrary organizations
535 // of 2d array data to be used.
536 //
537 // The subrange of indices used for contouring is kx to lx in the x
538 // direction and from ky to ly in the y direction. The array of contour
539 // levels is clevel(nlevel), and "pltr" is the name of a function which
540 // transforms array indicies into world coordinates.
541 //
542 // Note that the fortran-like minimum and maximum indices (kx, lx, ky, ly)
543 // are translated into more C-like ones. I've only kept them as they are
544 // for the plcontf() argument list because of backward compatibility.
545 //--------------------------------------------------------------------------
546 
547 void
548 plfcont( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ),
549  PLPointer f2eval_data,
550  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
551  PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel,
552  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
553  PLPointer pltr_data )
554 {
555  PLINT i, **ipts;
556 
557  if ( kx < 1 || kx >= lx )
558  {
559  plabort( "plfcont: indices must satisfy 1 <= kx <= lx <= nx" );
560  return;
561  }
562  if ( ky < 1 || ky >= ly )
563  {
564  plabort( "plfcont: indices must satisfy 1 <= ky <= ly <= ny" );
565  return;
566  }
567 
568  if ( ( ipts = (PLINT **) malloc( (size_t) nx * sizeof ( PLINT * ) ) ) == NULL )
569  {
570  plexit( "plfcont: Insufficient memory" );
571  }
572 
573  for ( i = 0; i < nx; i++ )
574  {
575  if ( ( ipts[i] = (PLINT *) malloc( (size_t) ny * sizeof ( PLINT * ) ) ) == NULL )
576  {
577  plexit( "plfcont: Insufficient memory" );
578  }
579  }
580 
581  for ( i = 0; i < nlevel; i++ )
582  {
583  plcntr( f2eval, f2eval_data,
584  nx, ny, kx - 1, lx - 1, ky - 1, ly - 1, clevel[i], ipts,
585  pltr, pltr_data );
586 
587  if ( error )
588  {
589  error = 0;
590  goto done;
591  }
592  }
593 
594 done:
595  for ( i = 0; i < nx; i++ )
596  {
597  free( (void *) ipts[i] );
598  }
599  free( (void *) ipts );
600 }
601 
602 //--------------------------------------------------------------------------
603 // void plcntr()
604 //
605 // The contour for a given level is drawn here. Note iscan has nx
606 // elements. ixstor and iystor each have nstor elements.
607 //--------------------------------------------------------------------------
608 
609 static void
610 plcntr( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ),
611  PLPointer f2eval_data,
612  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
613  PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts,
614  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
615  PLPointer pltr_data )
616 {
617  PLINT kcol, krow, lastindex;
618  PLFLT distance;
619  PLFLT save_def, save_scale;
620 
621  char flabel[30];
622  plgchr( &save_def, &save_scale );
623  save_scale = save_scale / save_def;
624 
625  cont_new_store( flev );
626 
627  // format contour label for plptex and define the font height of the labels
628  plfloatlabel( flev, flabel, 30 );
629  plschr( 0.0, contlabel_size );
630 
631  // Clear array for traversed squares
632  for ( kcol = kx; kcol < lx; kcol++ )
633  {
634  for ( krow = ky; krow < ly; krow++ )
635  {
636  ipts[kcol][krow] = 0;
637  }
638  }
639 
640 
641  for ( krow = ky; krow < ly; krow++ )
642  {
643  for ( kcol = kx; kcol < lx; kcol++ )
644  {
645  if ( ipts[kcol][krow] == 0 )
646  {
647  // Follow and draw a contour
648  pldrawcn( f2eval, f2eval_data,
649  nx, ny, kx, lx, ky, ly, flev, flabel, kcol, krow,
650  0.0, 0.0, -2, ipts, &distance, &lastindex,
651  pltr, pltr_data );
652 
653  if ( error )
654  return;
655  }
656  }
657  }
658  plschr( save_def, save_scale );
659 }
660 
661 //--------------------------------------------------------------------------
662 // void pldrawcn()
663 //
664 // Follow and draw a contour.
665 //--------------------------------------------------------------------------
666 
667 static void
668 pldrawcn( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ),
669  PLPointer f2eval_data,
670  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
671  PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow,
672  PLFLT lastx, PLFLT lasty, PLINT startedge, PLINT **ipts,
673  PLFLT *distance, PLINT *lastindex,
674  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
675  PLPointer pltr_data )
676 {
677  PLFLT f[4];
678  PLFLT px[4], py[4], locx[4], locy[4];
679  PLINT iedge[4];
680  PLINT i, j, k, num, first, inext, kcolnext, krownext, sfi, sfj;
681 
682 
683  ( *pltr )( kcol, krow + 1, &px[0], &py[0], pltr_data );
684  ( *pltr )( kcol, krow, &px[1], &py[1], pltr_data );
685  ( *pltr )( kcol + 1, krow, &px[2], &py[2], pltr_data );
686  ( *pltr )( kcol + 1, krow + 1, &px[3], &py[3], pltr_data );
687 
688  f[0] = f2eval( kcol, krow + 1, f2eval_data ) - flev;
689  f[1] = f2eval( kcol, krow, f2eval_data ) - flev;
690  f[2] = f2eval( kcol + 1, krow, f2eval_data ) - flev;
691  f[3] = f2eval( kcol + 1, krow + 1, f2eval_data ) - flev;
692 
693  for ( i = 0, j = 1; i < 4; i++, j = ( j + 1 ) % 4 )
694  {
695  // Use intermediates to avoid possible floating point
696  // under / over flow during multiplication.
697  sfi = ( f[i] > 0.0 ) ? 1 : ( ( f[i] < 0.0 ) ? -1 : 0 );
698  sfj = ( f[j] > 0.0 ) ? 1 : ( ( f[j] < 0.0 ) ? -1 : 0 );
699  iedge[i] = ( sfi * sfj > 0 ) ? -1 : ( ( sfi * sfj < 0 ) ? 1 : 0 );
700  }
701 
702  // Mark this square as done
703  ipts[kcol][krow] = 1;
704 
705  // Check if no contour has been crossed i.e. iedge[i] = -1
706  if ( ( iedge[0] == -1 ) && ( iedge[1] == -1 ) && ( iedge[2] == -1 )
707  && ( iedge[3] == -1 ) )
708  return;
709 
710  // Check if this is a completely flat square - in which case
711  // ignore it
712  if ( ( f[0] == 0.0 ) && ( f[1] == 0.0 ) && ( f[2] == 0.0 ) &&
713  ( f[3] == 0.0 ) )
714  return;
715 
716  // Calculate intersection points
717  num = 0;
718  if ( startedge < 0 )
719  {
720  first = 1;
721  }
722  else
723  {
724  locx[num] = lastx;
725  locy[num] = lasty;
726  num++;
727  first = 0;
728  }
729  for ( k = 0, i = ( startedge < 0 ? 0 : startedge ); k < 4; k++, i = ( i + 1 ) % 4 )
730  {
731  if ( i == startedge )
732  continue;
733 
734  // If the contour is an edge check it hasn't already been done
735  if ( f[i] == 0.0 && f[( i + 1 ) % 4] == 0.0 )
736  {
737  kcolnext = kcol;
738  krownext = krow;
739  if ( i == 0 )
740  kcolnext--;
741  if ( i == 1 )
742  krownext--;
743  if ( i == 2 )
744  kcolnext++;
745  if ( i == 3 )
746  krownext++;
747  if ( ( kcolnext < kx ) || ( kcolnext >= lx ) ||
748  ( krownext < ky ) || ( krownext >= ly ) ||
749  ( ipts[kcolnext][krownext] == 1 ) )
750  continue;
751  }
752  if ( ( iedge[i] == 1 ) || ( f[i] == 0.0 ) )
753  {
754  j = ( i + 1 ) % 4;
755  if ( f[i] != 0.0 )
756  {
757  locx[num] = ( px[i] * fabs( f[j] ) + px[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] );
758  locy[num] = ( py[i] * fabs( f[j] ) + py[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] );
759  }
760  else
761  {
762  locx[num] = px[i];
763  locy[num] = py[i];
764  }
765  // If this is the start of the contour then move to the point
766  if ( first == 1 )
767  {
768  cont_mv_store( locx[num], locy[num] );
769  first = 0;
770  *distance = 0;
771  *lastindex = 0;
772  }
773  else
774  {
775  // Link to the next point on the contour
776  if ( contlabel_active )
777  pl_drawcontlabel( locx[num], locy[num], flabel, distance, lastindex );
778  else
779  cont_xy_store( locx[num], locy[num] );
780  // Need to follow contour into next grid box
781  // Easy case where contour does not pass through corner
782  if ( f[i] != 0.0 )
783  {
784  kcolnext = kcol;
785  krownext = krow;
786  inext = ( i + 2 ) % 4;
787  if ( i == 0 )
788  kcolnext--;
789  if ( i == 1 )
790  krownext--;
791  if ( i == 2 )
792  kcolnext++;
793  if ( i == 3 )
794  krownext++;
795  if ( ( kcolnext >= kx ) && ( kcolnext < lx ) &&
796  ( krownext >= ky ) && ( krownext < ly ) &&
797  ( ipts[kcolnext][krownext] == 0 ) )
798  {
799  pldrawcn( f2eval, f2eval_data,
800  nx, ny, kx, lx, ky, ly, flev, flabel,
801  kcolnext, krownext,
802  locx[num], locy[num], inext, ipts,
803  distance, lastindex,
804  pltr, pltr_data );
805  }
806  }
807  // Hard case where contour passes through corner
808  // This is still not perfect - it may lose the contour
809  // which won't upset the contour itself (we can find it
810  // again later) but might upset the labelling
811  else
812  {
813  kcolnext = kcol;
814  krownext = krow;
815  inext = ( i + 2 ) % 4;
816  if ( i == 0 )
817  {
818  kcolnext--; krownext++;
819  }
820  if ( i == 1 )
821  {
822  krownext--; kcolnext--;
823  }
824  if ( i == 2 )
825  {
826  kcolnext++; krownext--;
827  }
828  if ( i == 3 )
829  {
830  krownext++; kcolnext++;
831  }
832  if ( ( kcolnext >= kx ) && ( kcolnext < lx ) &&
833  ( krownext >= ky ) && ( krownext < ly ) &&
834  ( ipts[kcolnext][krownext] == 0 ) )
835  {
836  pldrawcn( f2eval, f2eval_data,
837  nx, ny, kx, lx, ky, ly, flev, flabel,
838  kcolnext, krownext,
839  locx[num], locy[num], inext, ipts,
840  distance, lastindex,
841  pltr, pltr_data );
842  }
843  }
844  if ( first == 1 )
845  {
846  // Move back to first point
847  cont_mv_store( locx[num], locy[num] );
848  first = 0;
849  *distance = 0;
850  *lastindex = 0;
851  first = 0;
852  }
853  else
854  {
855  first = 1;
856  }
857  num++;
858  }
859  }
860  }
861 }
862 
863 //--------------------------------------------------------------------------
864 // pltr0()
865 //
866 // Identity transformation.
867 //--------------------------------------------------------------------------
868 
869 void
870 pltr0( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer PL_UNUSED( pltr_data ) )
871 {
872  *tx = x;
873  *ty = y;
874 }
875 
876 //--------------------------------------------------------------------------
877 // pltr1()
878 //
879 // Does linear interpolation from singly dimensioned coord arrays.
880 //
881 // Just abort for now if coordinates are out of bounds (don't think it's
882 // possible, but if so we could use linear extrapolation).
883 //--------------------------------------------------------------------------
884 
885 void
886 pltr1( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data )
887 {
888  PLINT ul, ur, vl, vr;
889  PLFLT du, dv;
890  PLFLT xl, xr, yl, yr;
891 
892  PLcGrid *grid = (PLcGrid *) pltr_data;
893  PLFLT *xg = grid->xg;
894  PLFLT *yg = grid->yg;
895  PLINT nx = grid->nx;
896  PLINT ny = grid->ny;
897 
898  ul = (PLINT) x;
899  ur = ul + 1;
900  du = x - ul;
901 
902  vl = (PLINT) y;
903  vr = vl + 1;
904  dv = y - vl;
905 
906  if ( x < 0 || x > nx - 1 || y < 0 || y > ny - 1 )
907  {
908  plexit( "pltr1: Invalid coordinates" );
909  }
910 
911 // Look up coordinates in row-dominant array.
912 // Have to handle right boundary specially -- if at the edge, we'd better
913 // not reference the out of bounds point.
914 //
915 
916  xl = xg[ul];
917  yl = yg[vl];
918 
919  if ( ur == nx )
920  {
921  *tx = xl;
922  }
923  else
924  {
925  xr = xg[ur];
926  *tx = xl * ( 1 - du ) + xr * du;
927  }
928  if ( vr == ny )
929  {
930  *ty = yl;
931  }
932  else
933  {
934  yr = yg[vr];
935  *ty = yl * ( 1 - dv ) + yr * dv;
936  }
937 }
938 
939 //--------------------------------------------------------------------------
940 // pltr2()
941 //
942 // Does linear interpolation from doubly dimensioned coord arrays (column
943 // dominant, as per normal C 2d arrays).
944 //
945 // This routine includes lots of checks for out of bounds. This would occur
946 // occasionally due to some bugs in the contour plotter (now fixed). If an
947 // out of bounds coordinate is obtained, the boundary value is provided
948 // along with a warning. These checks should stay since no harm is done if
949 // if everything works correctly.
950 //--------------------------------------------------------------------------
951 
952 void
953 pltr2( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data )
954 {
955  PLINT ul, ur, vl, vr;
956  PLFLT du, dv;
957  PLFLT xll, xlr, xrl, xrr;
958  PLFLT yll, ylr, yrl, yrr;
959  PLFLT xmin, xmax, ymin, ymax;
960 
961  PLcGrid2 *grid = (PLcGrid2 *) pltr_data;
962  PLFLT **xg = grid->xg;
963  PLFLT **yg = grid->yg;
964  PLINT nx = grid->nx;
965  PLINT ny = grid->ny;
966 
967  ul = (PLINT) x;
968  ur = ul + 1;
969  du = x - ul;
970 
971  vl = (PLINT) y;
972  vr = vl + 1;
973  dv = y - vl;
974 
975  xmin = 0;
976  xmax = nx - 1;
977  ymin = 0;
978  ymax = ny - 1;
979 
980  if ( x < xmin || x > xmax || y < ymin || y > ymax )
981  {
982  plwarn( "pltr2: Invalid coordinates" );
983  if ( x < xmin )
984  {
985  if ( y < ymin )
986  {
987  *tx = xg[0][0];
988  *ty = yg[0][0];
989  }
990  else if ( y > ymax )
991  {
992  *tx = xg[0][ny - 1];
993  *ty = yg[0][ny - 1];
994  }
995  else
996  {
997  xll = xg[0][vl];
998  yll = yg[0][vl];
999  xlr = xg[0][vr];
1000  ylr = yg[0][vr];
1001 
1002  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1003  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1004  }
1005  }
1006  else if ( x > xmax )
1007  {
1008  if ( y < ymin )
1009  {
1010  *tx = xg[nx - 1][0];
1011  *ty = yg[nx - 1][0];
1012  }
1013  else if ( y > ymax )
1014  {
1015  *tx = xg[nx - 1][ny - 1];
1016  *ty = yg[nx - 1][ny - 1];
1017  }
1018  else
1019  {
1020  xll = xg[nx - 1][vl];
1021  yll = yg[nx - 1][vl];
1022  xlr = xg[nx - 1][vr];
1023  ylr = yg[nx - 1][vr];
1024 
1025  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1026  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1027  }
1028  }
1029  else
1030  {
1031  if ( y < ymin )
1032  {
1033  xll = xg[ul][0];
1034  xrl = xg[ur][0];
1035  yll = yg[ul][0];
1036  yrl = yg[ur][0];
1037 
1038  *tx = xll * ( 1 - du ) + xrl * ( du );
1039  *ty = yll * ( 1 - du ) + yrl * ( du );
1040  }
1041  else if ( y > ymax )
1042  {
1043  xlr = xg[ul][ny - 1];
1044  xrr = xg[ur][ny - 1];
1045  ylr = yg[ul][ny - 1];
1046  yrr = yg[ur][ny - 1];
1047 
1048  *tx = xlr * ( 1 - du ) + xrr * ( du );
1049  *ty = ylr * ( 1 - du ) + yrr * ( du );
1050  }
1051  }
1052  }
1053 
1054 // Normal case.
1055 // Look up coordinates in row-dominant array.
1056 // Have to handle right boundary specially -- if at the edge, we'd
1057 // better not reference the out of bounds point.
1058 //
1059 
1060  else
1061  {
1062  xll = xg[ul][vl];
1063  yll = yg[ul][vl];
1064 
1065  // ur is out of bounds
1066 
1067  if ( ur == nx && vr < ny )
1068  {
1069  xlr = xg[ul][vr];
1070  ylr = yg[ul][vr];
1071 
1072  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1073  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1074  }
1075 
1076  // vr is out of bounds
1077 
1078  else if ( ur < nx && vr == ny )
1079  {
1080  xrl = xg[ur][vl];
1081  yrl = yg[ur][vl];
1082 
1083  *tx = xll * ( 1 - du ) + xrl * ( du );
1084  *ty = yll * ( 1 - du ) + yrl * ( du );
1085  }
1086 
1087  // both ur and vr are out of bounds
1088 
1089  else if ( ur == nx && vr == ny )
1090  {
1091  *tx = xll;
1092  *ty = yll;
1093  }
1094 
1095  // everything in bounds
1096 
1097  else
1098  {
1099  xrl = xg[ur][vl];
1100  xlr = xg[ul][vr];
1101  xrr = xg[ur][vr];
1102 
1103  yrl = yg[ur][vl];
1104  ylr = yg[ul][vr];
1105  yrr = yg[ur][vr];
1106 
1107  *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
1108  xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
1109 
1110  *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
1111  yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
1112  }
1113  }
1114 }
1115 
1116 //--------------------------------------------------------------------------
1117 // pltr2p()
1118 //
1119 // Just like pltr2() but uses pointer arithmetic to get coordinates from 2d
1120 // grid tables. This form of grid tables is compatible with those from
1121 // PLplot 4.0. The grid data must be pointed to by a PLcGrid structure.
1122 //--------------------------------------------------------------------------
1123 
1124 void
1125 pltr2p( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data )
1126 {
1127  PLINT ul, ur, vl, vr;
1128  PLFLT du, dv;
1129  PLFLT xll, xlr, xrl, xrr;
1130  PLFLT yll, ylr, yrl, yrr;
1131  PLFLT xmin, xmax, ymin, ymax;
1132 
1133  PLcGrid *grid = (PLcGrid *) pltr_data;
1134  PLFLT *xg = grid->xg;
1135  PLFLT *yg = grid->yg;
1136  PLINT nx = grid->nx;
1137  PLINT ny = grid->ny;
1138 
1139  ul = (PLINT) x;
1140  ur = ul + 1;
1141  du = x - ul;
1142 
1143  vl = (PLINT) y;
1144  vr = vl + 1;
1145  dv = y - vl;
1146 
1147  xmin = 0;
1148  xmax = nx - 1;
1149  ymin = 0;
1150  ymax = ny - 1;
1151 
1152  if ( x < xmin || x > xmax || y < ymin || y > ymax )
1153  {
1154  plwarn( "pltr2p: Invalid coordinates" );
1155  if ( x < xmin )
1156  {
1157  if ( y < ymin )
1158  {
1159  *tx = *xg;
1160  *ty = *yg;
1161  }
1162  else if ( y > ymax )
1163  {
1164  *tx = *( xg + ( ny - 1 ) );
1165  *ty = *( yg + ( ny - 1 ) );
1166  }
1167  else
1168  {
1169  ul = 0;
1170  xll = *( xg + ul * ny + vl );
1171  yll = *( yg + ul * ny + vl );
1172  xlr = *( xg + ul * ny + vr );
1173  ylr = *( yg + ul * ny + vr );
1174 
1175  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1176  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1177  }
1178  }
1179  else if ( x > xmax )
1180  {
1181  if ( y < ymin )
1182  {
1183  *tx = *( xg + ( ny - 1 ) * nx );
1184  *ty = *( yg + ( ny - 1 ) * nx );
1185  }
1186  else if ( y > ymax )
1187  {
1188  *tx = *( xg + ( ny - 1 ) + ( nx - 1 ) * ny );
1189  *ty = *( yg + ( ny - 1 ) + ( nx - 1 ) * ny );
1190  }
1191  else
1192  {
1193  ul = nx - 1;
1194  xll = *( xg + ul * ny + vl );
1195  yll = *( yg + ul * ny + vl );
1196  xlr = *( xg + ul * ny + vr );
1197  ylr = *( yg + ul * ny + vr );
1198 
1199  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1200  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1201  }
1202  }
1203  else
1204  {
1205  if ( y < ymin )
1206  {
1207  vl = 0;
1208  xll = *( xg + ul * ny + vl );
1209  xrl = *( xg + ur * ny + vl );
1210  yll = *( yg + ul * ny + vl );
1211  yrl = *( yg + ur * ny + vl );
1212 
1213  *tx = xll * ( 1 - du ) + xrl * ( du );
1214  *ty = yll * ( 1 - du ) + yrl * ( du );
1215  }
1216  else if ( y > ymax )
1217  {
1218  vr = ny - 1;
1219  xlr = *( xg + ul * ny + vr );
1220  xrr = *( xg + ur * ny + vr );
1221  ylr = *( yg + ul * ny + vr );
1222  yrr = *( yg + ur * ny + vr );
1223 
1224  *tx = xlr * ( 1 - du ) + xrr * ( du );
1225  *ty = ylr * ( 1 - du ) + yrr * ( du );
1226  }
1227  }
1228  }
1229 
1230 // Normal case.
1231 // Look up coordinates in row-dominant array.
1232 // Have to handle right boundary specially -- if at the edge, we'd better
1233 // not reference the out of bounds point.
1234 //
1235 
1236  else
1237  {
1238  xll = *( xg + ul * ny + vl );
1239  yll = *( yg + ul * ny + vl );
1240 
1241  // ur is out of bounds
1242 
1243  if ( ur == nx && vr < ny )
1244  {
1245  xlr = *( xg + ul * ny + vr );
1246  ylr = *( yg + ul * ny + vr );
1247 
1248  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1249  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1250  }
1251 
1252  // vr is out of bounds
1253 
1254  else if ( ur < nx && vr == ny )
1255  {
1256  xrl = *( xg + ur * ny + vl );
1257  yrl = *( yg + ur * ny + vl );
1258 
1259  *tx = xll * ( 1 - du ) + xrl * ( du );
1260  *ty = yll * ( 1 - du ) + yrl * ( du );
1261  }
1262 
1263  // both ur and vr are out of bounds
1264 
1265  else if ( ur == nx && vr == ny )
1266  {
1267  *tx = xll;
1268  *ty = yll;
1269  }
1270 
1271  // everything in bounds
1272 
1273  else
1274  {
1275  xrl = *( xg + ur * ny + vl );
1276  xlr = *( xg + ul * ny + vr );
1277  xrr = *( xg + ur * ny + vr );
1278 
1279  yrl = *( yg + ur * ny + vl );
1280  ylr = *( yg + ul * ny + vr );
1281  yrr = *( yg + ur * ny + vr );
1282 
1283  *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
1284  xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
1285 
1286  *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
1287  yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
1288  }
1289  }
1290 }
1291 
1292 //--------------------------------------------------------------------------
1293 // pltr2f()
1294 //
1295 // Does linear interpolation from doubly dimensioned coord arrays
1296 // (row dominant, i.e. Fortran ordering).
1297 //
1298 // This routine includes lots of checks for out of bounds. This would
1299 // occur occasionally due to a bug in the contour plotter that is now fixed.
1300 // If an out of bounds coordinate is obtained, the boundary value is provided
1301 // along with a warning. These checks should stay since no harm is done if
1302 // if everything works correctly.
1303 //--------------------------------------------------------------------------
1304 
1305 void
1306 pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data )
1307 {
1308  PLINT ul, ur, vl, vr;
1309  PLFLT du, dv;
1310  PLFLT xll, xlr, xrl, xrr;
1311  PLFLT yll, ylr, yrl, yrr;
1312  PLFLT xmin, xmax, ymin, ymax;
1313 
1314  PLcGrid *cgrid = (PLcGrid *) pltr_data;
1315  PLFLT *xg = cgrid->xg;
1316  PLFLT *yg = cgrid->yg;
1317  PLINT nx = cgrid->nx;
1318  PLINT ny = cgrid->ny;
1319 
1320  ul = (PLINT) x;
1321  ur = ul + 1;
1322  du = x - ul;
1323 
1324  vl = (PLINT) y;
1325  vr = vl + 1;
1326  dv = y - vl;
1327 
1328  xmin = 0;
1329  xmax = nx - 1;
1330  ymin = 0;
1331  ymax = ny - 1;
1332 
1333  if ( x < xmin || x > xmax || y < ymin || y > ymax )
1334  {
1335  plwarn( "pltr2f: Invalid coordinates" );
1336 
1337  if ( x < xmin )
1338  {
1339  if ( y < ymin )
1340  {
1341  *tx = *xg;
1342  *ty = *yg;
1343  }
1344  else if ( y > ymax )
1345  {
1346  *tx = *( xg + ( ny - 1 ) * nx );
1347  *ty = *( yg + ( ny - 1 ) * nx );
1348  }
1349  else
1350  {
1351  ul = 0;
1352  xll = *( xg + ul + vl * nx );
1353  yll = *( yg + ul + vl * nx );
1354  xlr = *( xg + ul + vr * nx );
1355  ylr = *( yg + ul + vr * nx );
1356 
1357  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1358  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1359  }
1360  }
1361  else if ( x > xmax )
1362  {
1363  if ( y < ymin )
1364  {
1365  *tx = *( xg + ( nx - 1 ) );
1366  *ty = *( yg + ( nx - 1 ) );
1367  }
1368  else if ( y > ymax )
1369  {
1370  *tx = *( xg + ( nx - 1 ) + ( ny - 1 ) * nx );
1371  *ty = *( yg + ( nx - 1 ) + ( ny - 1 ) * nx );
1372  }
1373  else
1374  {
1375  ul = nx - 1;
1376  xll = *( xg + ul + vl * nx );
1377  yll = *( yg + ul + vl * nx );
1378  xlr = *( xg + ul + vr * nx );
1379  ylr = *( yg + ul + vr * nx );
1380 
1381  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1382  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1383  }
1384  }
1385  else
1386  {
1387  if ( y < ymin )
1388  {
1389  vl = 0;
1390  xll = *( xg + ul + vl * nx );
1391  xrl = *( xg + ur + vl * nx );
1392  yll = *( yg + ul + vl * nx );
1393  yrl = *( yg + ur + vl * nx );
1394 
1395  *tx = xll * ( 1 - du ) + xrl * ( du );
1396  *ty = yll * ( 1 - du ) + yrl * ( du );
1397  }
1398  else if ( y > ymax )
1399  {
1400  vr = ny - 1;
1401  xlr = *( xg + ul + vr * nx );
1402  xrr = *( xg + ur + vr * nx );
1403  ylr = *( yg + ul + vr * nx );
1404  yrr = *( yg + ur + vr * nx );
1405 
1406  *tx = xlr * ( 1 - du ) + xrr * ( du );
1407  *ty = ylr * ( 1 - du ) + yrr * ( du );
1408  }
1409  }
1410  }
1411 
1412 // Normal case.
1413 // Look up coordinates in row-dominant array.
1414 // Have to handle right boundary specially -- if at the edge, we'd
1415 // better not reference the out of bounds point.
1416 
1417  else
1418  {
1419  xll = *( xg + ul + vl * nx );
1420  yll = *( yg + ul + vl * nx );
1421 
1422 // ur is out of bounds
1423 
1424  if ( ur == nx && vr < ny )
1425  {
1426  xlr = *( xg + ul + vr * nx );
1427  ylr = *( yg + ul + vr * nx );
1428 
1429  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1430  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1431  }
1432 
1433 // vr is out of bounds
1434 
1435  else if ( ur < nx && vr == ny )
1436  {
1437  xrl = *( xg + ur + vl * nx );
1438  yrl = *( yg + ur + vl * nx );
1439 
1440  *tx = xll * ( 1 - du ) + xrl * ( du );
1441  *ty = yll * ( 1 - du ) + yrl * ( du );
1442  }
1443 
1444 // both ur and vr are out of bounds
1445 
1446  else if ( ur == nx && vr == ny )
1447  {
1448  *tx = xll;
1449  *ty = yll;
1450  }
1451 
1452 // everything in bounds
1453 
1454  else
1455  {
1456  xrl = *( xg + ur + vl * nx );
1457  xlr = *( xg + ul + vr * nx );
1458  xrr = *( xg + ur + vr * nx );
1459 
1460  yrl = *( yg + ur + vl * nx );
1461  ylr = *( yg + ul + vr * nx );
1462  yrr = *( yg + ur + vr * nx );
1463  *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
1464  xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
1465 
1466  *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
1467  yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
1468  }
1469  }
1470 }