Menu
▾
▴
[Plplot-cvs] SF.net SVN: plplot:[10864] trunk
|
From: <ai...@us...> - 2010-03-13 23:37:36
|
Revision: 10864
http://plplot.svn.sourceforge.net/plplot/?rev=10864&view=rev
Author: airwin
Date: 2010-03-13 23:37:28 +0000 (Sat, 13 Mar 2010)
Log Message:
-----------
Apply patch from David MacMahon to support arbitrary storage of 2D
user data.
This is a roll-up of all previous PLplot patches related to supporting
arbitrary storage of 2D user data. This patch is based on (and should
apply cleanly to) svn trunk r10859.
Adds support for arbitrary storage of 2D user data. This is very
similar to the technique employed by some existing functions (e.g.
plfcont and plfshade) that use "evaluator" functions to access 2D user
data that is stored in an arbtrary format. The new approach extends
the concept of a user-supplied (or predefined) "evaluator" function to a
group of user-supplied (or predefined) "operator" functions. The
operator functions provide for various operations on the arbitrarily
stored 2D data including: get, set, +=, -=, *=, /=, isnan, minmax, and
f2eval.
To facilitate the passing of an entire family of operator functions (via
function pointers), a plf2ops_t structure is defined to contain a
pointer to each type of operator function. Predefined operator
functions are defined for several common 2D data storage techniques.
Variables (of type plf2ops_t) containing function pointers for these
operator functions are also defined.
New variants of functions that accept 2D data are created. The new
variants accept the 2D data as two parameters: a pointer to a plf2ops_t
structure containing (pointers to) suitable operator functions and a
PLPointer to the actual 2D data store. Existing functions that accept
2D data are modified to simply pass their parameters to the
corresponding new variant of the function, along with a pointer to the
suitable predefined plf2ops_t stucture of operator function pointers.
The list of functions for which new variants are created is:
c_plimage, c_plimagefr, c_plmesh, c_plmeshc, c_plot3d, c_plot3dc,
c_plot3dcl, c_plshade1, c_plshades, c_plsurf3d, and c_plsurf3dl, and
c_plgriddata. The new variants are named the same as their
corresponding existing function except that the "c_" prefix is changed
to "plf" (e.g. the new variant of c_plmesh is called plfmesh).
Adds plfvect declaration to plplot.h and changes the names (and only the
names) of some plfvect arguments to make them slightly clearer. In
order to maintain backwards API compatibility, this function and the
other existing functions that use "evaluator" functions are NOT changed
to use the new operator functions.
Makes plplot.h and libplplot consistent vis-a-vis pltr0f and pltr2d.
Moves the definitions of pltr2f (already declared in plplot.h) from the
sccont.c files of the FORTRAN 77 and Fortran 95 bindings into plcont.c.
Removes pltr0f declaration from plplot.h.
Changes x08c.c to demonstrate use of new support for arbitrary storage
of 2D data arrays. Shows how to do surface plots with the following
four types of 2D data arrays:
1) PLFLT z[nx][ny];
2) PLfGrid2 z;
3) PLFLT z[nx*ny]; /* row major order */
4) PLFLT z[nx*ny]; /* column major order */
---
bindings/f77/sccont.c | 182 ---------------------
bindings/f95/sccont.c | 182 ---------------------
examples/c/x08c.c | 45 ++++--
include/plplot.h | 238 ++++++++++++++++++++++++++--
src/CMakeLists.txt | 1 +
src/plcont.c | 207 +++++++++++++++++++++++-
src/plf2ops.c | 426 +++++++++++++++++++++++++++++++++++++++++++++++++
src/plgridd.c | 120 ++++++++------
src/plimage.c | 49 +++++--
src/plot3d.c | 212 ++++++++++++++++---------
src/plshade.c | 81 ++++++++--
src/plvect.c | 17 +--
12 files changed, 1200 insertions(+), 560 deletions(-)
Modified Paths:
--------------
trunk/bindings/f77/sccont.c
trunk/bindings/f95/sccont.c
trunk/examples/c/x08c.c
trunk/include/plplot.h
trunk/src/CMakeLists.txt
trunk/src/plcont.c
trunk/src/plfill.c
trunk/src/plgridd.c
trunk/src/plimage.c
trunk/src/plot3d.c
trunk/src/plshade.c
trunk/src/plvect.c
Modified: trunk/bindings/f77/sccont.c
===================================================================
--- trunk/bindings/f77/sccont.c 2010-03-13 19:39:09 UTC (rev 10863)
+++ trunk/bindings/f77/sccont.c 2010-03-13 23:37:28 UTC (rev 10864)
@@ -40,188 +40,6 @@
}
/*----------------------------------------------------------------------*\
- * pltr2f()
- *
- * Does linear interpolation from doubly dimensioned coord arrays
- * (row dominant, i.e. Fortran ordering).
- *
- * This routine includes lots of checks for out of bounds. This would
- * occur occasionally due to a bug in the contour plotter that is now fixed.
- * If an out of bounds coordinate is obtained, the boundary value is provided
- * along with a warning. These checks should stay since no harm is done if
- * if everything works correctly.
- \*----------------------------------------------------------------------*/
-
-void
-pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data )
-{
- PLINT ul, ur, vl, vr;
- PLFLT du, dv;
- PLFLT xll, xlr, xrl, xrr;
- PLFLT yll, ylr, yrl, yrr;
- PLFLT xmin, xmax, ymin, ymax;
-
- PLcGrid *cgrid = (PLcGrid *) pltr_data;
- PLFLT *xg = cgrid->xg;
- PLFLT *yg = cgrid->yg;
- PLINT nx = cgrid->nx;
- PLINT ny = cgrid->ny;
-
- ul = (PLINT) x;
- ur = ul + 1;
- du = x - ul;
-
- vl = (PLINT) y;
- vr = vl + 1;
- dv = y - vl;
-
- xmin = 0;
- xmax = nx - 1;
- ymin = 0;
- ymax = ny - 1;
-
- if ( x < xmin || x > xmax || y < ymin || y > ymax )
- {
- plwarn( "pltr2f: Invalid coordinates" );
-
- if ( x < xmin )
- {
- if ( y < ymin )
- {
- *tx = *xg;
- *ty = *yg;
- }
- else if ( y > ymax )
- {
- *tx = *( xg + ( ny - 1 ) * nx );
- *ty = *( yg + ( ny - 1 ) * nx );
- }
- else
- {
- ul = 0;
- xll = *( xg + ul + vl * nx );
- yll = *( yg + ul + vl * nx );
- xlr = *( xg + ul + vr * nx );
- ylr = *( yg + ul + vr * nx );
-
- *tx = xll * ( 1 - dv ) + xlr * ( dv );
- *ty = yll * ( 1 - dv ) + ylr * ( dv );
- }
- }
- else if ( x > xmax )
- {
- if ( y < ymin )
- {
- *tx = *( xg + ( nx - 1 ) );
- *ty = *( yg + ( nx - 1 ) );
- }
- else if ( y > ymax )
- {
- *tx = *( xg + ( nx - 1 ) + ( ny - 1 ) * nx );
- *ty = *( yg + ( nx - 1 ) + ( ny - 1 ) * nx );
- }
- else
- {
- ul = nx - 1;
- xll = *( xg + ul + vl * nx );
- yll = *( yg + ul + vl * nx );
- xlr = *( xg + ul + vr * nx );
- ylr = *( yg + ul + vr * nx );
-
- *tx = xll * ( 1 - dv ) + xlr * ( dv );
- *ty = yll * ( 1 - dv ) + ylr * ( dv );
- }
- }
- else
- {
- if ( y < ymin )
- {
- vl = 0;
- xll = *( xg + ul + vl * nx );
- xrl = *( xg + ur + vl * nx );
- yll = *( yg + ul + vl * nx );
- yrl = *( yg + ur + vl * nx );
-
- *tx = xll * ( 1 - du ) + xrl * ( du );
- *ty = yll * ( 1 - du ) + yrl * ( du );
- }
- else if ( y > ymax )
- {
- vr = ny - 1;
- xlr = *( xg + ul + vr * nx );
- xrr = *( xg + ur + vr * nx );
- ylr = *( yg + ul + vr * nx );
- yrr = *( yg + ur + vr * nx );
-
- *tx = xlr * ( 1 - du ) + xrr * ( du );
- *ty = ylr * ( 1 - du ) + yrr * ( du );
- }
- }
- }
-
-/* Normal case.
- * Look up coordinates in row-dominant array.
- * Have to handle right boundary specially -- if at the edge, we'd
- * better not reference the out of bounds point. */
-
- else
- {
- xll = *( xg + ul + vl * nx );
- yll = *( yg + ul + vl * nx );
-
-/* ur is out of bounds */
-
- if ( ur == nx && vr < ny )
- {
- xlr = *( xg + ul + vr * nx );
- ylr = *( yg + ul + vr * nx );
-
- *tx = xll * ( 1 - dv ) + xlr * ( dv );
- *ty = yll * ( 1 - dv ) + ylr * ( dv );
- }
-
-/* vr is out of bounds */
-
- else if ( ur < nx && vr == ny )
- {
- xrl = *( xg + ur + vl * nx );
- yrl = *( yg + ur + vl * nx );
-
- *tx = xll * ( 1 - du ) + xrl * ( du );
- *ty = yll * ( 1 - du ) + yrl * ( du );
- }
-
-/* both ur and vr are out of bounds */
-
- else if ( ur == nx && vr == ny )
- {
- *tx = xll;
- *ty = yll;
- }
-
-/* everything in bounds */
-
- else
- {
- xrl = *( xg + ur + vl * nx );
- xlr = *( xg + ul + vr * nx );
- xrr = *( xg + ur + vr * nx );
-
- yrl = *( yg + ur + vl * nx );
- ylr = *( yg + ul + vr * nx );
- yrr = *( yg + ur + vr * nx );
-/* INDENT OFF */
- *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
- xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
-
- *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
- yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
-/* INDENT ON */
- }
- }
-}
-
-/*----------------------------------------------------------------------*\
* Contour plotter front-ends.
* These specify the row-dominant function evaluator in the plfcont
* argument list. NO TRANSPOSE IS NECESSARY. The routines are as follows:
Modified: trunk/bindings/f95/sccont.c
===================================================================
--- trunk/bindings/f95/sccont.c 2010-03-13 19:39:09 UTC (rev 10863)
+++ trunk/bindings/f95/sccont.c 2010-03-13 23:37:28 UTC (rev 10864)
@@ -40,188 +40,6 @@
}
/*----------------------------------------------------------------------*\
- * pltr2f()
- *
- * Does linear interpolation from doubly dimensioned coord arrays
- * (row dominant, i.e. Fortran ordering).
- *
- * This routine includes lots of checks for out of bounds. This would
- * occur occasionally due to a bug in the contour plotter that is now fixed.
- * If an out of bounds coordinate is obtained, the boundary value is provided
- * along with a warning. These checks should stay since no harm is done if
- * if everything works correctly.
- \*----------------------------------------------------------------------*/
-
-void
-pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data )
-{
- PLINT ul, ur, vl, vr;
- PLFLT du, dv;
- PLFLT xll, xlr, xrl, xrr;
- PLFLT yll, ylr, yrl, yrr;
- PLFLT xmin, xmax, ymin, ymax;
-
- PLcGrid *cgrid = (PLcGrid *) pltr_data;
- PLFLT *xg = cgrid->xg;
- PLFLT *yg = cgrid->yg;
- PLINT nx = cgrid->nx;
- PLINT ny = cgrid->ny;
-
- ul = (PLINT) x;
- ur = ul + 1;
- du = x - ul;
-
- vl = (PLINT) y;
- vr = vl + 1;
- dv = y - vl;
-
- xmin = 0;
- xmax = nx - 1;
- ymin = 0;
- ymax = ny - 1;
-
- if ( x < xmin || x > xmax || y < ymin || y > ymax )
- {
- plwarn( "pltr2f: Invalid coordinates" );
-
- if ( x < xmin )
- {
- if ( y < ymin )
- {
- *tx = *xg;
- *ty = *yg;
- }
- else if ( y > ymax )
- {
- *tx = *( xg + ( ny - 1 ) * nx );
- *ty = *( yg + ( ny - 1 ) * nx );
- }
- else
- {
- ul = 0;
- xll = *( xg + ul + vl * nx );
- yll = *( yg + ul + vl * nx );
- xlr = *( xg + ul + vr * nx );
- ylr = *( yg + ul + vr * nx );
-
- *tx = xll * ( 1 - dv ) + xlr * ( dv );
- *ty = yll * ( 1 - dv ) + ylr * ( dv );
- }
- }
- else if ( x > xmax )
- {
- if ( y < ymin )
- {
- *tx = *( xg + ( nx - 1 ) );
- *ty = *( yg + ( nx - 1 ) );
- }
- else if ( y > ymax )
- {
- *tx = *( xg + ( nx - 1 ) + ( ny - 1 ) * nx );
- *ty = *( yg + ( nx - 1 ) + ( ny - 1 ) * nx );
- }
- else
- {
- ul = nx - 1;
- xll = *( xg + ul + vl * nx );
- yll = *( yg + ul + vl * nx );
- xlr = *( xg + ul + vr * nx );
- ylr = *( yg + ul + vr * nx );
-
- *tx = xll * ( 1 - dv ) + xlr * ( dv );
- *ty = yll * ( 1 - dv ) + ylr * ( dv );
- }
- }
- else
- {
- if ( y < ymin )
- {
- vl = 0;
- xll = *( xg + ul + vl * nx );
- xrl = *( xg + ur + vl * nx );
- yll = *( yg + ul + vl * nx );
- yrl = *( yg + ur + vl * nx );
-
- *tx = xll * ( 1 - du ) + xrl * ( du );
- *ty = yll * ( 1 - du ) + yrl * ( du );
- }
- else if ( y > ymax )
- {
- vr = ny - 1;
- xlr = *( xg + ul + vr * nx );
- xrr = *( xg + ur + vr * nx );
- ylr = *( yg + ul + vr * nx );
- yrr = *( yg + ur + vr * nx );
-
- *tx = xlr * ( 1 - du ) + xrr * ( du );
- *ty = ylr * ( 1 - du ) + yrr * ( du );
- }
- }
- }
-
-/* Normal case.
- * Look up coordinates in row-dominant array.
- * Have to handle right boundary specially -- if at the edge, we'd
- * better not reference the out of bounds point. */
-
- else
- {
- xll = *( xg + ul + vl * nx );
- yll = *( yg + ul + vl * nx );
-
-/* ur is out of bounds */
-
- if ( ur == nx && vr < ny )
- {
- xlr = *( xg + ul + vr * nx );
- ylr = *( yg + ul + vr * nx );
-
- *tx = xll * ( 1 - dv ) + xlr * ( dv );
- *ty = yll * ( 1 - dv ) + ylr * ( dv );
- }
-
-/* vr is out of bounds */
-
- else if ( ur < nx && vr == ny )
- {
- xrl = *( xg + ur + vl * nx );
- yrl = *( yg + ur + vl * nx );
-
- *tx = xll * ( 1 - du ) + xrl * ( du );
- *ty = yll * ( 1 - du ) + yrl * ( du );
- }
-
-/* both ur and vr are out of bounds */
-
- else if ( ur == nx && vr == ny )
- {
- *tx = xll;
- *ty = yll;
- }
-
-/* everything in bounds */
-
- else
- {
- xrl = *( xg + ur + vl * nx );
- xlr = *( xg + ul + vr * nx );
- xrr = *( xg + ur + vr * nx );
-
- yrl = *( yg + ur + vl * nx );
- ylr = *( yg + ul + vr * nx );
- yrr = *( yg + ur + vr * nx );
-/* INDENT OFF */
- *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
- xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
-
- *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
- yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
-/* INDENT ON */
- }
- }
-}
-
-/*----------------------------------------------------------------------*\
* Contour plotter front-ends.
* These specify the row-dominant function evaluator in the plfcont
* argument list. NO TRANSPOSE IS NECESSARY. The routines are as follows:
Modified: trunk/examples/c/x08c.c
===================================================================
--- trunk/examples/c/x08c.c 2010-03-13 19:39:09 UTC (rev 10863)
+++ trunk/examples/c/x08c.c 2010-03-13 23:37:28 UTC (rev 10864)
@@ -26,6 +26,10 @@
#include "plcdemos.h"
+/* plexit not declared in public header! */
+PLDLLIMPEXP void
+plexit( const char *errormsg );
+
#define XPTS 35 /* Data points in x */
#define YPTS 46 /* Data points in y */
@@ -120,14 +124,15 @@
int
main( int argc, const char *argv[] )
{
- int i, j, k;
- PLFLT *x, *y, **z;
- PLFLT xx, yy, r;
- PLINT ifshade;
- PLFLT zmin, zmax, step;
- PLFLT clevel[LEVELS];
- PLINT nlevel = LEVELS;
- int rosen = 1;
+ int i, j, k;
+ PLFLT *x, *y, **z, *z_row_major, *z_col_major;
+ PLfGrid2 grid_c, grid_row_major, grid_col_major;
+ PLFLT xx, yy, r;
+ PLINT ifshade;
+ PLFLT zmin, zmax, step;
+ PLFLT clevel[LEVELS];
+ PLINT nlevel = LEVELS;
+ int rosen = 1;
/* Parse and process command line arguments */
plMergeOpts( options, "x08c options", NULL );
@@ -145,7 +150,17 @@
y = (PLFLT *) calloc( YPTS, sizeof ( PLFLT ) );
plAlloc2dGrid( &z, XPTS, YPTS );
+ z_row_major = (PLFLT *) malloc( XPTS * YPTS * sizeof ( PLFLT ) );
+ z_col_major = (PLFLT *) malloc( XPTS * YPTS * sizeof ( PLFLT ) );
+ if ( !z_row_major || !z_col_major )
+ plexit( "Memory allocation error" );
+ grid_c.f = z;
+ grid_row_major.f = (PLFLT **) z_row_major;
+ grid_col_major.f = (PLFLT **) z_col_major;
+ grid_c.nx = grid_row_major.nx = grid_col_major.nx = XPTS;
+ grid_c.ny = grid_row_major.ny = grid_col_major.ny = YPTS;
+
for ( i = 0; i < XPTS; i++ )
{
x[i] = ( (double) ( i - ( XPTS / 2 ) ) / (double) ( XPTS / 2 ) );
@@ -169,6 +184,7 @@
if ( rosen )
{
z[i][j] = pow( 1. - xx, 2. ) + 100. * pow( yy - pow( xx, 2. ), 2. );
+
/* The log argument may be zero for just the right grid. */
if ( z[i][j] > 0. )
z[i][j] = log( z[i][j] );
@@ -180,6 +196,9 @@
r = sqrt( xx * xx + yy * yy );
z[i][j] = exp( -r * r ) * cos( 2.0 * M_PI * r );
}
+
+ z_row_major[i * YPTS + j] = z[i][j];
+ z_col_major[i + XPTS * j] = z[i][j];
}
}
@@ -213,22 +232,22 @@
if ( ifshade == 0 ) /* diffuse light surface plot */
{
cmap1_init( 1 );
- plsurf3d( x, y, z, XPTS, YPTS, 0, NULL, 0 );
+ plfsurf3d( x, y, plf2ops_c(), (PLPointer) z, XPTS, YPTS, 0, NULL, 0 );
}
else if ( ifshade == 1 ) /* magnitude colored plot */
{
cmap1_init( 0 );
- plsurf3d( x, y, z, XPTS, YPTS, MAG_COLOR, NULL, 0 );
+ plfsurf3d( x, y, plf2ops_grid_c(), ( PLPointer ) & grid_c, XPTS, YPTS, MAG_COLOR, NULL, 0 );
}
else if ( ifshade == 2 ) /* magnitude colored plot with faceted squares */
{
cmap1_init( 0 );
- plsurf3d( x, y, z, XPTS, YPTS, MAG_COLOR | FACETED, NULL, 0 );
+ plfsurf3d( x, y, plf2ops_grid_row_major(), ( PLPointer ) & grid_row_major, XPTS, YPTS, MAG_COLOR | FACETED, NULL, 0 );
}
else /* magnitude colored plot with contours */
{
cmap1_init( 0 );
- plsurf3d( x, y, z, XPTS, YPTS, MAG_COLOR | SURF_CONT | BASE_CONT, clevel, nlevel );
+ plfsurf3d( x, y, plf2ops_grid_col_major(), ( PLPointer ) & grid_col_major, XPTS, YPTS, MAG_COLOR | SURF_CONT | BASE_CONT, clevel, nlevel );
}
}
}
@@ -238,6 +257,8 @@
free( (void *) x );
free( (void *) y );
plFree2dGrid( z, XPTS, YPTS );
+ free( (void *) z_row_major );
+ free( (void *) z_col_major );
plend();
Modified: trunk/include/plplot.h
===================================================================
--- trunk/include/plplot.h 2010-03-13 19:39:09 UTC (rev 10863)
+++ trunk/include/plplot.h 2010-03-13 23:37:28 UTC (rev 10864)
@@ -191,7 +191,6 @@
typedef PLINT PLBOOL;
/* For passing user data, as with X's XtPointer */
-
typedef void* PLPointer;
/*--------------------------------------------------------------------------*\
@@ -527,6 +526,42 @@
PLINT nValues;
} PLAttribute;
+/*
+ * typedefs for access methods for arbitrary (i.e. user defined) data storage
+ */
+
+/*
+ * This type of struct holds pointers to functions that are used to get, set,
+ * modify, and test individual 2-D data points referenced by a PLPointer. How
+ * the PLPointer is used depends entirely on the functions that implement the
+ * various operations. Certain common data representations have predefined
+ * instances of this structure prepopulated with pointers to predefined
+ * functions.
+ */
+
+typedef struct
+{
+ PLFLT ( *get )( PLPointer p, PLINT ix, PLINT iy );
+ PLFLT ( *set )( PLPointer p, PLINT ix, PLINT iy, PLFLT z );
+ PLFLT ( *add )( PLPointer p, PLINT ix, PLINT iy, PLFLT z );
+ PLFLT ( *sub )( PLPointer p, PLINT ix, PLINT iy, PLFLT z );
+ PLFLT ( *mul )( PLPointer p, PLINT ix, PLINT iy, PLFLT z );
+ PLFLT ( *div )( PLPointer p, PLINT ix, PLINT iy, PLFLT z );
+ PLINT ( *is_nan )( PLPointer p, PLINT ix, PLINT iy );
+ void ( *minmax )( PLPointer p, PLINT nx, PLINT ny, PLFLT *zmim, PLFLT *zmax );
+ /*
+ * f2eval is backwards compatible signature for "f2eval" functions that
+ * existed before plf2ops "operator function families" were used.
+ */
+ PLFLT ( *f2eval )( PLINT ix, PLINT iy, PLPointer p );
+} plf2ops_t;
+
+/*
+ * A typedef to facilitate declaration of a pointer to a plfops_t structure.
+ */
+
+typedef plf2ops_t * PLF2OPS;
+
/*--------------------------------------------------------------------------*\
* BRAINDEAD-ness
*
@@ -1079,6 +1114,11 @@
PLFLT *xg, PLINT nptsx, PLFLT *yg, PLINT nptsy,
PLFLT **zg, PLINT type, PLFLT data );
+PLDLLIMPEXP void
+plfgriddata( PLFLT *x, PLFLT *y, PLFLT *z, PLINT npts,
+ PLFLT *xg, PLINT nptsx, PLFLT *yg, PLINT nptsy,
+ PLF2OPS zops, PLPointer zgp, PLINT type, PLFLT data );
+
/* type of gridding algorithm for plgriddata() */
#define GRID_CSA 1 /* Bivariate Cubic Spline approximation */
@@ -1199,12 +1239,24 @@
PLDLLIMPEXP void
c_plmesh( PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt );
+/* Like plmesh, but uses an evaluator function to access z data from zp */
+
+PLDLLIMPEXP void
+plfmesh( PLFLT *x, PLFLT *y, PLF2OPS zops, PLPointer zp,
+ PLINT nx, PLINT ny, PLINT opt );
+
/* Plots a mesh representation of the function z[x][y] with contour */
PLDLLIMPEXP void
c_plmeshc( PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt,
PLFLT *clevel, PLINT nlevel );
+/* Like plmeshc, but uses an evaluator function to access z data from zp */
+
+PLDLLIMPEXP void
+plfmeshc( PLFLT *x, PLFLT *y, PLF2OPS zops, PLPointer zp,
+ PLINT nx, PLINT ny, PLINT opt, PLFLT *clevel, PLINT nlevel );
+
/* Creates a new stream and makes it the default. */
PLDLLIMPEXP void
@@ -1228,6 +1280,12 @@
c_plot3d( PLFLT *x, PLFLT *y, PLFLT **z,
PLINT nx, PLINT ny, PLINT opt, PLBOOL side );
+/* Like plot3d, but uses an evaluator function to access z data from zp */
+
+PLDLLIMPEXP void
+plfplot3d( PLFLT *x, PLFLT *y, PLF2OPS zops, PLPointer zp,
+ PLINT nx, PLINT ny, PLINT opt, PLBOOL side );
+
/* Plots a 3-d representation of the function z[x][y] with contour. */
PLDLLIMPEXP void
@@ -1235,6 +1293,12 @@
PLINT nx, PLINT ny, PLINT opt,
PLFLT *clevel, PLINT nlevel );
+/* Like plot3dc, but uses an evaluator function to access z data from zp */
+
+PLDLLIMPEXP void
+plfplot3dc( PLFLT *x, PLFLT *y, PLF2OPS zops, PLPointer zp,
+ PLINT nx, PLINT ny, PLINT opt, PLFLT *clevel, PLINT nlevel );
+
/* Plots a 3-d representation of the function z[x][y] with contour and
* y index limits. */
@@ -1244,6 +1308,14 @@
PLFLT *clevel, PLINT nlevel,
PLINT ixstart, PLINT ixn, PLINT *indexymin, PLINT*indexymax );
+/* Like plot3dcl, but uses an evaluator function to access z data from zp */
+
+PLDLLIMPEXP void
+plfplot3dcl( PLFLT *x, PLFLT *y, PLF2OPS zops, PLPointer zp,
+ PLINT nx, PLINT ny, PLINT opt,
+ PLFLT *clevel, PLINT nlevel,
+ PLINT ixstart, PLINT ixn, PLINT *indexymin, PLINT *indexymax );
+
/*
* definitions for the opt argument in plot3dc() and plsurf3d()
*
@@ -1505,6 +1577,16 @@
PLPointer pltr_data );
PLDLLIMPEXP void
+plfshades( PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
+ PLINT ( *defined )( PLFLT, PLFLT ),
+ PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
+ PLFLT *clevel, PLINT nlevel, PLINT fill_width,
+ PLINT cont_color, PLINT cont_width,
+ void ( *fill )( PLINT, PLFLT *, PLFLT * ), PLINT rectangular,
+ void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
+ PLPointer pltr_data );
+
+PLDLLIMPEXP void
plfshade( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ),
PLPointer f2eval_data,
PLFLT ( *c2eval )( PLINT, PLINT, PLPointer ),
@@ -1519,6 +1601,18 @@
void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
PLPointer pltr_data );
+PLDLLIMPEXP void
+plfshade1( PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
+ PLINT ( *defined )( PLFLT, PLFLT ),
+ PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
+ PLFLT shade_min, PLFLT shade_max,
+ PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
+ PLINT min_color, PLINT min_width,
+ PLINT max_color, PLINT max_width,
+ void ( *fill )( PLINT, PLFLT *, PLFLT * ), PLINT rectangular,
+ void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
+ PLPointer pltr_data );
+
/* Setup a user-provided custom labeling function */
PLDLLIMPEXP void
@@ -1621,6 +1715,18 @@
void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
PLPointer pltr_data );
+/*
+ * Like plimagefr, but uses an evaluator function to access image data from
+ * idatap. getminmax is only used if zmin == zmax.
+ */
+
+PLDLLIMPEXP void
+plfimagefr( PLF2OPS idataops, PLPointer idatap, PLINT nx, PLINT ny,
+ PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax, PLFLT zmin, PLFLT zmax,
+ PLFLT valuemin, PLFLT valuemax,
+ void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
+ PLPointer pltr_data );
+
/* plots a 2d image (or a matrix too large for plshade() ) - colors
* automatically scaled */
@@ -1629,6 +1735,16 @@
PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax, PLFLT zmin, PLFLT zmax,
PLFLT Dxmin, PLFLT Dxmax, PLFLT Dymin, PLFLT Dymax );
+/*
+ * Like plimage, but uses an operator functions to access image data from
+ * idatap.
+ */
+
+PLDLLIMPEXP void
+plfimage( PLF2OPS idataops, PLPointer idatap, PLINT nx, PLINT ny,
+ PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax, PLFLT zmin, PLFLT zmax,
+ PLFLT Dxmin, PLFLT Dxmax, PLFLT Dymin, PLFLT Dymax );
+
/* Set up a new line style */
PLDLLIMPEXP void
@@ -1640,6 +1756,12 @@
c_plsurf3d( PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny,
PLINT opt, PLFLT *clevel, PLINT nlevel );
+/* Like plsurf3d, but uses an evaluator function to access z data from zp */
+
+PLDLLIMPEXP void
+plfsurf3d( PLFLT *x, PLFLT *y, PLF2OPS zops, PLPointer zp,
+ PLINT nx, PLINT ny, PLINT opt, PLFLT *clevel, PLINT nlevel );
+
/* Plots the 3d surface representation of the function z[x][y] with y
* index limits. */
@@ -1648,7 +1770,14 @@
PLINT opt, PLFLT *clevel, PLINT nlevel,
PLINT ixstart, PLINT ixn, PLINT *indexymin, PLINT*indexymax );
+/* Like plsurf3dl, but uses an evaluator function to access z data from zp */
+
PLDLLIMPEXP void
+plfsurf3dl( PLFLT *x, PLFLT *y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
+ PLINT opt, PLFLT *clevel, PLINT nlevel,
+ PLINT ixstart, PLINT ixn, PLINT *indexymin, PLINT *indexymax );
+
+PLDLLIMPEXP void
c_plsvect( PLFLT *arrowx, PLFLT *arrowy, PLINT npts, PLBOOL fill );
/* Sets the edges of the viewport to the specified absolute coordinates */
@@ -1707,7 +1836,18 @@
void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
PLPointer pltr_data );
+/*
+ * Routine to plot a vector array with arbitrary coordinate
+ * and vector transformations
+ */
PLDLLIMPEXP void
+plfvect( PLFLT ( *getuv )( PLINT, PLINT, PLPointer ),
+ PLPointer up, PLPointer vp,
+ PLINT nx, PLINT ny, PLFLT scale,
+ void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
+ PLPointer pltr_data );
+
+PLDLLIMPEXP void
c_plvpas( PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax, PLFLT aspect );
/* Creates a viewport with the specified normalized subpage coordinates. */
@@ -1826,33 +1966,105 @@
PLDLLIMPEXP void
pltr2p( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data );
-/* Identity transformation for plots from Fortran. */
-
-PLDLLIMPEXP void
-pltr0f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data );
-
/* Does linear interpolation from doubly dimensioned coord arrays */
/* (row dominant, i.e. Fortran ordering). */
PLDLLIMPEXP void
pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data );
-/* Function evaluators */
+/*
+ * Returns a pointer to a plf2ops_t stucture with pointers to functions for
+ * accessing 2-D data referenced as (PLFLT **), such as the C variable z
+ * declared as...
+ *
+ * PLFLT z[nx][ny];
+ */
-/* Does a lookup from a 2d function array. Array is of type (PLFLT **), */
-/* and is column dominant (normal C ordering). */
+PLDLLIMPEXP PLF2OPS
+plf2ops_c();
+/*
+ * Returns a pointer to a plf2ops_t stucture with pointers to functions for accessing 2-D data
+ * referenced as (PLfGrid2 *), where the PLfGrid2's "f" is treated as type
+ * (PLFLT **).
+ */
+
+PLDLLIMPEXP PLF2OPS
+plf2ops_grid_c();
+
+/*
+ * Returns a pointer to a plf2ops_t stucture with pointers to functions for
+ * accessing 2-D data stored in (PLfGrid2 *), with the PLfGrid2's "f" field
+ * treated as type (PLFLT *) pointing to 2-D data stored in row-major order.
+ * In the context of plotting, it might be easier to think of it as "X-major"
+ * order. In this ordering, values for a single X index are stored in
+ * consecutive memory locations.
+ */
+
+PLDLLIMPEXP PLF2OPS
+plf2ops_grid_row_major();
+
+/*
+ * Returns a pointer to a plf2ops_t stucture with pointers to functions for
+ * accessing 2-D data stored in (PLfGrid2 *), with the PLfGrid2's "f" field
+ * treated as type (PLFLT *) pointing to 2-D data stored in column-major order.
+ * In the context of plotting, it might be easier to think of it as "Y-major"
+ * order. In this ordering, values for a single Y index are stored in
+ * consecutive memory locations.
+ */
+
+PLDLLIMPEXP PLF2OPS
+plf2ops_grid_col_major();
+
+
+/* Function evaluators (Should these be deprecated in favor of plf2ops?) */
+
+/*
+ * Does a lookup from a 2d function array. plf2eval_data is treated as type
+ * (PLFLT **) and data for (ix,iy) is returned from...
+ *
+ * plf2eval_data[ix][iy];
+ */
+
PLDLLIMPEXP PLFLT
+plf2eval1( PLINT ix, PLINT iy, PLPointer plf2eval_data );
+
+/*
+ * Does a lookup from a 2d function array. plf2eval_data is treated as type
+ * (PLfGrid2 *) and data for (ix,iy) is returned from...
+ *
+ * plf2eval_data->f[ix][iy];
+ */
+
+PLDLLIMPEXP PLFLT
plf2eval2( PLINT ix, PLINT iy, PLPointer plf2eval_data );
-/* Does a lookup from a 2d function array. Array is of type (PLFLT *), */
-/* and is column dominant (normal C ordering). */
+/*
+ * Does a lookup from a 2d function array. plf2eval_data is treated as type
+ * (PLfGrid *) and data for (ix,iy) is returned from...
+ *
+ * plf2eval_data->f[ix * plf2eval_data->ny + iy];
+ *
+ * This is commonly called "row-major order", but in the context of plotting,
+ * it might be easier to think of it as "X-major order". In this ordering,
+ * values for a single X index are stored in consecutive memory locations.
+ * This is also known as C ordering.
+ */
PLDLLIMPEXP PLFLT
plf2eval( PLINT ix, PLINT iy, PLPointer plf2eval_data );
-/* Does a lookup from a 2d function array. Array is of type (PLFLT *), */
-/* and is row dominant (Fortran ordering). */
+/*
+ * Does a lookup from a 2d function array. plf2eval_data is treated as type
+ * (PLfGrid *) and data for (ix,iy) is returned from...
+ *
+ * plf2eval_data->f[ix + iy * plf2eval_data->nx];
+ *
+ * This is commonly called "column-major order", but in the context of
+ * plotting, it might be easier to think of it as "Y-major order". In this
+ * ordering, values for a single Y index are stored in consecutive memory
+ * locations. This is also known as FORTRAN ordering.
+ */
PLDLLIMPEXP PLFLT
plf2evalr( PLINT ix, PLINT iy, PLPointer plf2eval_data );
Modified: trunk/src/CMakeLists.txt
===================================================================
--- trunk/src/CMakeLists.txt 2010-03-13 19:39:09 UTC (rev 10863)
+++ trunk/src/CMakeLists.txt 2010-03-13 23:37:28 UTC (rev 10864)
@@ -31,6 +31,7 @@
plcvt.c
pldeprecated.c
pldtik.c
+ plf2ops.c
plfill.c
plfreetype.c
plgradient.c
Modified: trunk/src/plcont.c
===================================================================
--- trunk/src/plcont.c 2010-03-13 19:39:09 UTC (rev 10863)
+++ trunk/src/plcont.c 2010-03-13 23:37:28 UTC (rev 10864)
@@ -405,13 +405,29 @@
return ( ( y - plsc->wpyoff ) / plsc->wpyscl );
}
+/*--------------------------------------------------------------------------*\
+ * plf2eval1()
+ *
+ * Does a lookup from a 2d function array. Array is of type (PLFLT **),
+ * and is column dominant (normal C ordering).
+ \*--------------------------------------------------------------------------*/
+PLFLT
+plf2eval1( PLINT ix, PLINT iy, PLPointer plf2eval_data )
+{
+ PLFLT value;
+ PLFLT **z = (PLFLT **) plf2eval_data;
+ value = z[ix][iy];
+
+ return value;
+}
+
/*--------------------------------------------------------------------------*\
* plf2eval2()
*
- * Does a lookup from a 2d function array. Array is of type (PLFLT **),
- * and is column dominant (normal C ordering).
+ * Does a lookup from a 2d function array. plf2eval_data is treated as type
+ * (PLfGrid2 *).
\*--------------------------------------------------------------------------*/
PLFLT
@@ -501,8 +517,6 @@
void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
PLPointer pltr_data )
{
- PLfGrid2 grid;
-
if ( pltr == NULL )
{
/* If pltr is undefined, abort with an error. */
@@ -510,8 +524,7 @@
return;
}
- grid.f = f;
- plfcont( plf2eval2, ( PLPointer ) & grid,
+ plfcont( plf2eval1, ( PLPointer ) f,
nx, ny, kx, lx, ky, ly, clevel, nlevel,
pltr, pltr_data );
}
@@ -1261,3 +1274,185 @@
}
}
}
+
+/*----------------------------------------------------------------------*\
+ * pltr2f()
+ *
+ * Does linear interpolation from doubly dimensioned coord arrays
+ * (row dominant, i.e. Fortran ordering).
+ *
+ * This routine includes lots of checks for out of bounds. This would
+ * occur occasionally due to a bug in the contour plotter that is now fixed.
+ * If an out of bounds coordinate is obtained, the boundary value is provided
+ * along with a warning. These checks should stay since no harm is done if
+ * if everything works correctly.
+ \*----------------------------------------------------------------------*/
+
+void
+pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data )
+{
+ PLINT ul, ur, vl, vr;
+ PLFLT du, dv;
+ PLFLT xll, xlr, xrl, xrr;
+ PLFLT yll, ylr, yrl, yrr;
+ PLFLT xmin, xmax, ymin, ymax;
+
+ PLcGrid *cgrid = (PLcGrid *) pltr_data;
+ PLFLT *xg = cgrid->xg;
+ PLFLT *yg = cgrid->yg;
+ PLINT nx = cgrid->nx;
+ PLINT ny = cgrid->ny;
+
+ ul = (PLINT) x;
+ ur = ul + 1;
+ du = x - ul;
+
+ vl = (PLINT) y;
+ vr = vl + 1;
+ dv = y - vl;
+
+ xmin = 0;
+ xmax = nx - 1;
+ ymin = 0;
+ ymax = ny - 1;
+
+ if ( x < xmin || x > xmax || y < ymin || y > ymax )
+ {
+ plwarn( "pltr2f: Invalid coordinates" );
+
+ if ( x < xmin )
+ {
+ if ( y < ymin )
+ {
+ *tx = *xg;
+ *ty = *yg;
+ }
+ else if ( y > ymax )
+ {
+ *tx = *( xg + ( ny - 1 ) * nx );
+ *ty = *( yg + ( ny - 1 ) * nx );
+ }
+ else
+ {
+ ul = 0;
+ xll = *( xg + ul + vl * nx );
+ yll = *( yg + ul + vl * nx );
+ xlr = *( xg + ul + vr * nx );
+ ylr = *( yg + ul + vr * nx );
+
+ *tx = xll * ( 1 - dv ) + xlr * ( dv );
+ *ty = yll * ( 1 - dv ) + ylr * ( dv );
+ }
+ }
+ else if ( x > xmax )
+ {
+ if ( y < ymin )
+ {
+ *tx = *( xg + ( nx - 1 ) );
+ *ty = *( yg + ( nx - 1 ) );
+ }
+ else if ( y > ymax )
+ {
+ *tx = *( xg + ( nx - 1 ) + ( ny - 1 ) * nx );
+ *ty = *( yg + ( nx - 1 ) + ( ny - 1 ) * nx );
+ }
+ else
+ {
+ ul = nx - 1;
+ xll = *( xg + ul + vl * nx );
+ yll = *( yg + ul + vl * nx );
+ xlr = *( xg + ul + vr * nx );
+ ylr = *( yg + ul + vr * nx );
+
+ *tx = xll * ( 1 - dv ) + xlr * ( dv );
+ *ty = yll * ( 1 - dv ) + ylr * ( dv );
+ }
+ }
+ else
+ {
+ if ( y < ymin )
+ {
+ vl = 0;
+ xll = *( xg + ul + vl * nx );
+ xrl = *( xg + ur + vl * nx );
+ yll = *( yg + ul + vl * nx );
+ yrl = *( yg + ur + vl * nx );
+
+ *tx = xll * ( 1 - du ) + xrl * ( du );
+ *ty = yll * ( 1 - du ) + yrl * ( du );
+ }
+ else if ( y > ymax )
+ {
+ vr = ny - 1;
+ xlr = *( xg + ul + vr * nx );
+ xrr = *( xg + ur + vr * nx );
+ ylr = *( yg + ul + vr * nx );
+ yrr = *( yg + ur + vr * nx );
+
+ *tx = xlr * ( 1 - du ) + xrr * ( du );
+ *ty = ylr * ( 1 - du ) + yrr * ( du );
+ }
+ }
+ }
+
+/* Normal case.
+ * Look up coordinates in row-dominant array.
+ * Have to handle right boundary specially -- if at the edge, we'd
+ * better not reference the out of bounds point. */
+
+ else
+ {
+ xll = *( xg + ul + vl * nx );
+ yll = *( yg + ul + vl * nx );
+
+/* ur is out of bounds */
+
+ if ( ur == nx && vr < ny )
+ {
+ xlr = *( xg + ul + vr * nx );
+ ylr = *( yg + ul + vr * nx );
+
+ *tx = xll * ( 1 - dv ) + xlr * ( dv );
+ *ty = yll * ( 1 - dv ) + ylr * ( dv );
+ }
+
+/* vr is out of bounds */
+
+ else if ( ur < nx && vr == ny )
+ {
+ xrl = *( xg + ur + vl * nx );
+ yrl = *( yg + ur + vl * nx );
+
+ *tx = xll * ( 1 - du ) + xrl * ( du );
+ *ty = yll * ( 1 - du ) + yrl * ( du );
+ }
+
+/* both ur and vr are out of bounds */
+
+ else if ( ur == nx && vr == ny )
+ {
+ *tx = xll;
+ *ty = yll;
+ }
+
+/* everything in bounds */
+
+ else
+ {
+ xrl = *( xg + ur + vl * nx );
+ xlr = *( xg + ul + vr * nx );
+ xrr = *( xg + ur + vr * nx );
+
+ yrl = *( yg + ur + vl * nx );
+ ylr = *( yg + ul + vr * nx );
+ yrr = *( yg + ur + vr * nx );
+/* INDENT OFF */
+ *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
+ xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
+
+ *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
+ yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
+/* INDENT ON */
+ }
+ }
+}
Modified: trunk/src/plfill.c
===================================================================
--- trunk/src/plfill.c 2010-03-13 19:39:09 UTC (rev 10863)
+++ trunk/src/plfill.c 2010-03-13 23:37:28 UTC (rev 10864)
@@ -1895,9 +1895,10 @@
status = status | PL_NEAR_PARALLEL;
else
status = status | PL_PARALLEL;
- /* Choice of intersect is arbitrary in this case. Choose A1, A2,
- * B1, or B2 (in that order) if any of them lie inside or near
- * the other line segment. Otherwise, choose the average point. */
+ /* Choice of intersect is arbitrary in this case. Choose A1
+ * A2, B1, or B2 (in that order) if any of them lie inside or
+ * near the other line segment. Otherwise, choose the average
+ * point. */
if ( ( BETW_NBCC( xA1, xB1, xB2 ) && BETW_NBCC( yA1, yB1, yB2 ) ) )
{
fxintersect = xA1;
@@ -1945,11 +1946,11 @@
* Find out which. */
if ( fabs( fxintersect - xA1 ) <= PL_NBCC && fabs( fyintersect - yA1 ) <= PL_NBCC )
status = status | PL_NEAR_A1;
- else if ( fabs( fxintersect - xA2 ) <= PL_NBCC && fabs( fyintersect - yA2 ) <= PL_NBCC )
+ if ( fabs( fxintersect - xA2 ) <= PL_NBCC && fabs( fyintersect - yA2 ) <= PL_NBCC )
status = status | PL_NEAR_A2;
- else if ( fabs( fxintersect - xB1 ) <= PL_NBCC && fabs( fyintersect - yB1 ) <= PL_NBCC )
+ if ( fabs( fxintersect - xB1 ) <= PL_NBCC && fabs( fyintersect - yB1 ) <= PL_NBCC )
status = status | PL_NEAR_B1;
- else if ( fabs( fxintersect - xB2 ) <= PL_NBCC && fabs( fyintersect - yB2 ) <= PL_NBCC )
+ if ( fabs( fxintersect - xB2 ) <= PL_NBCC && fabs( fyintersect - yB2 ) <= PL_NBCC )
status = status | PL_NEAR_B2;
/* N.B. if none of the above conditions hold then status remains at
* zero to signal we have a definite crossing. */
@@ -2021,6 +2022,13 @@
int i1m1, i2, i2m1, ifnotcrossed;
int ifxsort, ifascend, count_crossings = 0, status = 0;
PLINT xintersect, yintersect;
+ PLINT ifnear_a1, ifnear_a2, ifnear_a, ifnear_b1, ifnear_b2, ifnear_b;
+ PLINT crossing_index[4];
+ /* The following are floating point to avoid integer overflows while
+ * still having exact arithmetic (except for very large values beyond
+ * the 53 bits of the significand). */
+ PLFLT xline, yline, xpoint0, ypoint0, xpoint1, ypoint1,
+ inprod0, inprod1, product_NBCC;
i1m1 = i1 - 1;
if ( i1m1 < 0 )
@@ -2053,7 +2061,124 @@
ifnotcrossed = notcrossed( &xintersect, &yintersect,
x1[i1m1], y1[i1m1], x1[i1], y1[i1],
x2[i2m1], y2[i2m1], x2[i2], y2[i2] );
+ if ( ifnotcrossed && !( ifnotcrossed & ( PL_PARALLEL | PL_NEAR_PARALLEL | PL_NOT_CROSSED ) ) )
+ {
+ /* intersection near at least one of the vertices, but not
+ * near parallel and not parallel */
+ ifnear_a1 = ifnotcrossed & PL_NEAR_A1;
+ ifnear_a2 = ifnotcrossed & PL_NEAR_A2;
+ ifnear_a = ifnotcrossed & ( PL_NEAR_A1 | PL_NEAR_A2 );
+ ifnear_b1 = ifnotcrossed & PL_NEAR_B1;
+ ifnear_b2 = ifnotcrossed & PL_NEAR_B2;
+ ifnear_b = ifnotcrossed & ( PL_NEAR_B1 | PL_NEAR_B2 );
+ if ( ( ifnear_a1 && ifnear_a2 ) || ( ifnear_b1 && ifnear_b2 ) ||
+ !( ifnear_a || ifnear_b ) )
+ plwarn( "number_crossings; internal logic error" );
+ if ( ifnear_a1 )
+ {
+ crossing_index[0] = i1;
+ crossing_index[2] = i1m1 - 1;
+ if ( crossing_index[2] < 0 )
+ crossing_index[2] += n1;
+ }
+ else if ( ifnear_a2 )
+ {
+ crossing_index[0] = i1m1;
+ crossing_index[2] = i1 + 1;
+ if ( crossing_index[2] >= n1 )
+ crossing_index[2] -= n1;
+ }
+ if ( ifnear_b1 )
+ {
+ crossing_index[1] = i2;
+ crossing_index[3] = i2m1 - 1;
+ if ( crossing_index[3] < 0 )
+ crossing_index[3] += n2;
+ }
+ else if ( ifnear_b2 )
+ {
+ crossing_index[1] = i2m1;
+ crossing_index[3] = i2 + 1;
+ if ( crossing_index[3] >= n2 )
+ crossing_index[3] -= n2;
+ }
+ /* Never count ifnear_a1 since that should be a duplicate
+ * of ifnear_a2 for the next i1 index. */
+ if ( ifnear_a2 && !ifnear_b )
+ {
+ /* Are the points defined by crossing_index[0] (point
+ * 0) and crossing_index[2] (point 1) on opposite
+ * sides of the line defined by the polygon 2 edge
+ * from i2m1 to i2? If so, indicate definite
+ * crossing. */
+ /* Use point i2m1 as origin of coordinate system and take
+ * dot products of both point 0 and point 1 vectors with
+ * perpendicular to point i2 vector. */
+ xline = -( y2[i2] - y2[i2m1] );
+ yline = x2[i2] - x2[i2m1];
+ xpoint0 = x1[crossing_index[0]] - x2[i2m1];
+ ypoint0 = y1[crossing_index[0]] - y2[i2m1];
+ xpoint1 = x1[crossing_index[2]] - x2[i2m1];
+ ypoint1 = y1[crossing_index[2]] - y2[i2m1];
+ inprod0 = xpoint0 * xline + ypoint0 * yline;
+ inprod1 = xpoint1 * xline + ypoint1 * yline;
+ /* maximum error of product of inner products for
+ * an error of PL_NBCC of the right sign for
+ * each original position. */
+ product_NBCC = 2. * PL_NBCC * (
+ ( fabs( xpoint0 ) + fabs( xline ) +
+ fabs( ypoint0 ) + fabs( yline ) ) * fabs( inprod1 ) +
+ ( fabs( xpoint1 ) + fabs( xline ) +
+ fabs( ypoint1 ) + fabs( yline ) ) * fabs( inprod0 ) );
+
+ /* Not crossed if sine 0 proportional to inprod0 has
+ * same sign (subject to error) as sine 1 proportional
+ * to inprod1. */
+ ifnotcrossed = inprod0 * inprod1 >= -product_NBCC;
+ }
+ /* Never count ifnear_b1 since that should be a duplicate
+ * of ifnear_b2 for the next i2 index. */
+ else if ( ifnear_b2 && !ifnear_a )
+ {
+ /* Are the points defined by crossing_index[1] (point
+ * 0) and crossing_index[3] (point 1) on opposite
+ * sides of the line defined by the polygon 1 edge
+ * from i1m1 to i1? If so, indicate definite
+ * crossing. */
+
+ /* Use point i1m1 as origin of coordinate system and take
+ * dot products of both point 0 and point 1 vectors with
+ * perpendicular to point i1 vector. */
+ xline = -( y1[i1] - y1[i1m1] );
+ yline = x1[i1] - x1[i1m1];
+ xpoint0 = x2[crossing_index[1]] - x1[i1m1];
+ ypoint0 = y2[crossing_index[1]] - y1[i1m1];
+ xpoint1 = x2[crossing_index[3]] - x1[i1m1];
+ ypoint1 = y2[crossing_index[3]] - y1[i1m1];
+ inprod0 = xpoint0 * xline + ypoint0 * yline;
+ inprod1 = xpoint1 * xline + ypoint1 * yline;
+ /* maximum error of product of inner products for
+ * an error of PL_NBCC of the right sign for
+ * each original position. */
+ product_NBCC = 2. * PL_NBCC * (
+ ( fabs( xpoint0 ) + fabs( xline ) +
+ fabs( ypoint0 ) + fabs( yline ) ) * fabs( inprod1 ) +
+ ( fabs( xpoint1 ) + fabs( xline ) +
+ fabs( ypoint1 ) + fabs( yline ) ) * fabs( inprod0 ) );
+
+ /* Not crossed if sine 0 proportional to inprod0 has
+ * same sign (subject to error) as sine 1 proportional
+ * to inprod1. */
+ ifnotcrossed = inprod0 * inprod1 >= -product_NBCC;
+ }
+ else if ( ifnear_a && ifnear_b )
+ {
+ plwarn( "notcrossing: coincident vertices not implemented yet." );
+ ifnotcrossed = 0;
+ }
+ }
+
if ( !ifnotcrossed )
{
count_crossings++;
Modified: trunk/src/plgridd.c
===================================================================
--- trunk/src/plgridd.c 2010-03-13 19:39:09 UTC (rev 10863)
+++ trunk/src/plgridd.c 2010-03-13 23:37:28 UTC (rev 10864)
@@ -38,33 +38,36 @@
/* forward declarations */
static void
grid_nnaidw( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg );
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp );
static void
grid_nnli( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg,
- PLFLT threshold );
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp, PLFLT threshold );
static void
grid_nnidw( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg,
- int knn_order );
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp, int knn_order );
#ifdef WITH_CSA
static void
grid_csa( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg );
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp );
#endif
#ifdef HAVE_QHULL
static void
grid_nni( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg,
- PLFLT wmin );
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp, PLFLT wmin );
static void
grid_dtli( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg );
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp );
#endif
static void
@@ -112,6 +115,14 @@
PLFLT *xg, PLINT nptsx, PLFLT *yg, PLINT nptsy,
PLFLT **zg, PLINT type, PLFLT data )
{
+ plfgriddata( x, y, z, npts, xg, nptsx, yg, nptsy, plf2ops_c(), (PLPointer) zg, type, data );
+}
+
+void
+plfgriddata( PLFLT *x, PLFLT *y, PLFLT *z, PLINT npts,
+ PLFLT *xg, PLINT nptsx, PLFLT *yg, PLINT nptsy,
+ PLF2OPS zops, PLPointer zgp, PLINT type, PLFLT data )
+{
int i, j;
if ( npts < 1 || nptsx < 1 || nptsy < 1 )
@@ -142,46 +153,46 @@
/* clear array to return */
for ( i = 0; i < nptsx; i++ )
for ( j = 0; j < nptsy; j++ )
- zg[i][j] = 0.; /* NaN signals a not processed grid point */
+ zops->set( zgp, i, j, 0.0 ); /* NaN signals a not processed grid point */
switch ( type )
{
case ( GRID_CSA ): /* Bivariate Cubic Spline Approximation */
#ifdef WITH_CSA
- grid_csa( x, y, z, npts, xg, nptsx, yg, nptsy, zg );
+ grid_csa( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
#else
plwarn( "plgriddata(): PLplot was configured to not use GRID_CSA.\n Reverting to GRID_NNAIDW." );
- grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zg );
+ grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
#endif
break;
case ( GRID_NNIDW ): /* Nearest Neighbors Inverse Distance Weighted */
- grid_nnidw( x, y, z, npts, xg, nptsx, yg, nptsy, zg, (int) data );
+ grid_nnidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, (int) data );
break;
case ( GRID_NNLI ): /* Nearest Neighbors Linear Interpolation */
- grid_nnli( x, y, z, npts, xg, nptsx, yg, nptsy, zg, data );
+ grid_nnli( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, data );
break;
case ( GRID_NNAIDW ): /* Nearest Neighbors "Around" Inverse Distance Weighted */
- grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zg );
+ grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
break;
case ( GRID_DTLI ): /* Delaunay Triangulation Linear Interpolation */
#ifdef HAVE_QHULL
- grid_dtli( x, y, z, npts, xg, nptsx, yg, nptsy, zg );
+ grid_dtli( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
#else
plwarn( "plgriddata(): you must have the Qhull library installed to use GRID_DTLI.\n Reverting to GRID_NNAIDW." );
- grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zg );
+ grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
#endif
break;
case ( GRID_NNI ): /* Natural Neighbors */
#ifdef HAVE_QHULL
- grid_nni( x, y, z, npts, xg, nptsx, yg, nptsy, zg, data );
+ grid_nni( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, data );
#else
plwarn( "plgriddata(): you must have the Qhull library installed to use GRID_NNI.\n Reverting to GRID_NNAIDW." );
- grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zg );
+ grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
#endif
break;
@@ -199,7 +210,8 @@
static void
grid_csa( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg )
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp )
{
PLFLT *xt, *yt, *zt;
point *pin, *pgrid, *pt;
@@ -248,8 +260,8 @@
{
for ( j = 0; j < nptsy; j++ )
{
- pt = &pgrid[j * nptsx + i];
- zg[i][j] = (PLFLT) pt->z;
+ pt = &pgrid[j * nptsx + i];
+ zops->set( zgp, i, j, (PLFLT) pt->z );
}
}
@@ -269,8 +281,8 @@
static void
grid_nnidw( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg,
- int knn_order )
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp, int knn_order )
{
int i, j, k;
PLFLT wi, nt;
@@ -300,8 +312,8 @@
if ( items[k].dist > md )
md = items[k].dist;
#endif
- zg[i][j] = 0.;
- nt = 0.;
+ zops->set( zgp, i, j, 0.0 );
+ nt = 0.;
for ( k = 0; k < knn_order; k++ )
{
@@ -313,13 +325,13 @@
#else
wi = 1. / ( items[k].dist * items[k].dist );
#endif
- zg[i][j] += wi * z[items[k].item];
- nt += wi;
+ zops->add( zgp, i, j, wi * z[items[k].item] );
+ nt += wi;
}
if ( nt != 0. )
- zg[i][j] /= nt;
+ zops->div( zgp, i, j, nt );
else
- zg[i][j] = NaN;
+ zops->set( zgp, i, j, NaN );
}
}
}
@@ -332,8 +344,8 @@
static void
grid_nnli( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg,
- PLFLT threshold )
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy,
+ PLF2OPS zops, PLPointer zgp, PLFLT threshold )
{
PLFLT xx[4], yy[4], zz[4], t, A, B, C, D, d1, d2, d3, max_thick;
int i, j, ii, excl, cnt, excl_item;
@@ -369,7 +381,7 @@
if ( d1 == 0. || d2 == 0. || d3 == 0. ) /* coincident points */
{
- zg[i][j] = NaN;
+ zops->set( zgp, i, j, NaN );
continue;
}
@@ -387,7 +399,7 @@
if ( ( d1 + d2 ) / d3 < threshold ) /* thin triangle! */
{
- zg[i][j] = NaN; /* deal with it latter */
+ zops->set( zgp, i, j, NaN ); /* deal with it later */
}
else /* calculate the plane passing through the three points */
@@ -398,7 +410,7 @@
D = -A * xx[0] - B * yy[0] - C * zz[0];
/* and interpolate (or extrapolate...) */
- zg[i][j] = -xg[i] * A / C - yg[j] * B / C - D / C;
+ zops->set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
}
}
}
@@ -417,7 +429,7 @@
{
for ( j = 0; j < nptsy; j++ )
{
- if ( isnan( zg[i][j] ) )
+ if ( zops->is_nan( zgp, i, j ) )
{
dist1( xg[i], yg[j], x, y, npts, 4 );
@@ -494,7 +506,7 @@
D = -A * xx[0] - B * yy[0] - C * zz[0];
/* and interpolate (or extrapolate...) */
- zg[i][j] = -xg[i] * A / C - yg[j] * B / C - D / C;
+ zops->set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
}
}
}
@@ -510,7 +522,7 @@
static void
grid_nnaidw( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg )
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLF2OPS zops, PLPointer zgp )
{
PLFLT d, nt;
int i, j, k;
@@ -520,21 +532,21 @@
for ( j = 0; j < nptsy; j++ )
{
dist2( xg[i], yg[j], x, y, npts );
- zg[i][j] = 0.;
- nt = 0.;
+ zops->set( zgp, i, j, 0. );
+ nt = 0.;
for ( k = 0; k < 4; k++ )
{
if ( items[k].item != -1 ) /* was found */
{
- d = 1. / ( items[k].dist * items[k].dist ); /* 1/square distance */
- zg[i][j] += d * z[items[k].item];
- nt += d;
+ d = 1. / ( items[k].dist * items[k].dist ); /* 1/square distance */
+ zops->add( zgp, i, j, d * z[items[k].item] );
+ nt += d;
}
}
if ( nt == 0. ) /* no points found?! */
- zg[i][j] = NaN;
+ zops->set( zgp, i, j, NaN );
else
- zg[i][j] /= nt;
+ zops->div( zgp, i, j, nt );
}
}
}
@@ -553,7 +565,7 @@
static void
grid_dtli( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg )
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLF2OPS zops, PLPointer zgp )
{
point *pin, *pgrid, *pt;
PLFLT *xt, *yt, *zt;
@@ -604,8 +616,8 @@
{
for ( j = 0; j < nptsy; j++ )
{
- pt = &pgrid[j * nptsx + i];
- zg[i][j] = (PLFLT) pt->z;
+ pt = &pgrid[j * nptsx + i];
+ zops->set( zgp, i, j, (PLFLT) pt->z );
}
}
@@ -622,7 +634,7 @@
static void
grid_nni( PLFLT *x, PLFLT *y, PLFLT *z, int npts,
- PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg,
+ PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLF2OPS zops, PLPointer zgp,
PLFLT wmin )
{
PLFLT *xt, *yt, *zt;
@@ -681,8 +693,8 @@
{
for ( j = 0; j < nptsy; j++ )
{
- pt = &pgrid[j * nptsx + i];
- zg[i][j] = (PLFLT) pt->z;
+ pt = &pgrid[j * nptsx + i];
+ zops->set( ...
[truncated message content] |
