35#ifdef HAS_LIBQHULL_INCLUDE
36#include <libqhull/qhull_a.h>
38#include <qhull/qhull_a.h>
83#define KNN_MAX_ORDER 100
123 plfgriddata( x, y, z, npts, xg, nptsx, yg, nptsy,
plf2ops_c(), (
PLPointer)
zg, type, data );
133 if ( npts < 1 || nptsx < 1 || nptsy < 1 )
135 plabort(
"plgriddata: Bad array dimensions" );
141 for ( i = 0; i < nptsx - 1; i++ )
143 if ( xg[i] >= xg[i + 1] )
145 plabort(
"plgriddata: xg array must be strictly increasing" );
149 for ( i = 0; i < nptsy - 1; i++ )
151 if ( yg[i] >= yg[i + 1] )
153 plabort(
"plgriddata: yg array must be strictly increasing" );
159 for ( i = 0; i < nptsx; i++ )
160 for ( j = 0; j < nptsy; j++ )
161 zops->set( zgp, i, j, 0.0 );
168 grid_csa( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
170 plwarn(
"plgriddata(): PLplot was configured to not use GRID_CSA.\n Reverting to GRID_NNAIDW." );
171 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
176 grid_nnidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, (
int) data );
180 grid_nnli( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, data );
184 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
189 grid_dtli( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
191 plwarn(
"plgriddata(): you must have the Qhull library installed to use GRID_DTLI.\n Reverting to GRID_NNAIDW." );
192 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
198 grid_nni( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, data );
200 plwarn(
"plgriddata(): you must have the Qhull library installed to use GRID_NNI.\n Reverting to GRID_NNAIDW." );
201 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
206 plabort(
"plgriddata: unknown algorithm type" );
227 if ( ( pin = (
point *) malloc( (
size_t) npts *
sizeof (
point ) ) ) == NULL )
229 plexit(
"grid_csa: Insufficient memory" );
236 for ( i = 0; i < npts; i++ )
238 pt->x = (double) *xt++;
239 pt->y = (double) *yt++;
240 pt->z = (double) *zt++;
244 nptsg = nptsx * nptsy;
245 if ( ( pgrid = (
point *) malloc( (
size_t) nptsg *
sizeof (
point ) ) ) == NULL )
247 plexit(
"grid_csa: Insufficient memory" );
252 for ( j = 0; j < nptsy; j++ )
255 for ( i = 0; i < nptsx; i++ )
257 pt->x = (double) *xt++;
258 pt->y = (double) *yt;
269 for ( i = 0; i < nptsx; i++ )
271 for ( j = 0; j < nptsy; j++ )
273 pt = &pgrid[j * nptsx + i];
274 zops->set( zgp, i, j, (
PLFLT)
pt->z );
302 plabort(
"plgriddata(): GRID_NNIDW: knn_order too big" );
306 if ( knn_order == 0 )
308 plwarn(
"plgriddata(): GRID_NNIDW: knn_order must be specified with 'data' arg. Using 15" );
312 for ( i = 0; i < nptsx; i++ )
314 for ( j = 0; j < nptsy; j++ )
316 dist1( xg[i], yg[j], x, y, npts, knn_order );
321 for ( k = 1; k < knn_order; k++ )
322 if (
items[k].dist > md )
325 zops->set( zgp, i, j, 0.0 );
328 for ( k = 0; k < knn_order; k++ )
330 if (
items[k].item == -1 )
333 wi = ( md -
items[k].dist ) / ( md *
items[k].dist );
338 zops->add( zgp, i, j, wi * z[
items[k].item] );
342 zops->div( zgp, i, j, nt );
344 zops->set( zgp, i, j,
NaN );
360 PLFLT xx[4], yy[4], zz[4], t, A, B, C, D, d1, d2, d3, max_thick;
361 int i, j, ii, excl, cnt, excl_item;
363 if ( threshold == 0. )
365 plwarn(
"plgriddata(): GRID_NNLI: threshold must be specified with 'data' arg. Using 1.001" );
368 else if ( threshold > 2. || threshold < 1. )
370 plabort(
"plgriddata(): GRID_NNLI: 1. < threshold < 2." );
374 for ( i = 0; i < nptsx; i++ )
376 for ( j = 0; j < nptsy; j++ )
378 dist1( xg[i], yg[j], x, y, npts, 3 );
381 for ( ii = 0; ii < 3; ii++ )
383 xx[ii] = x[
items[ii].item];
384 yy[ii] = y[
items[ii].item];
385 zz[ii] = z[
items[ii].item];
388 d1 = sqrt( ( xx[1] - xx[0] ) * ( xx[1] - xx[0] ) + ( yy[1] - yy[0] ) * ( yy[1] - yy[0] ) );
389 d2 = sqrt( ( xx[2] - xx[1] ) * ( xx[2] - xx[1] ) + ( yy[2] - yy[1] ) * ( yy[2] - yy[1] ) );
390 d3 = sqrt( ( xx[0] - xx[2] ) * ( xx[0] - xx[2] ) + ( yy[0] - yy[2] ) * ( yy[0] - yy[2] ) );
392 if ( d1 == 0. || d2 == 0. || d3 == 0. )
394 zops->set( zgp, i, j,
NaN );
401 t = d1; d1 = d2; d2 = t;
407 t = d2; d2 = d3; d3 = t;
410 if ( ( d1 + d2 ) / d3 < threshold )
412 zops->set( zgp, i, j,
NaN );
417 A = yy[0] * ( zz[1] - zz[2] ) + yy[1] * ( zz[2] - zz[0] ) + yy[2] * ( zz[0] - zz[1] );
418 B = zz[0] * ( xx[1] - xx[2] ) + zz[1] * ( xx[2] - xx[0] ) + zz[2] * ( xx[0] - xx[1] );
419 C = xx[0] * ( yy[1] - yy[2] ) + xx[1] * ( yy[2] - yy[0] ) + xx[2] * ( yy[0] - yy[1] );
420 D = -A * xx[0] - B * yy[0] - C * zz[0];
423 zops->set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
438 for ( i = 0; i < nptsx; i++ )
440 for ( j = 0; j < nptsy; j++ )
442 if ( zops->is_nan( zgp, i, j ) )
444 dist1( xg[i], yg[j], x, y, npts, 4 );
458 max_thick = 0.; excl_item = -1;
459 for ( excl = 0; excl < 4; excl++ )
463 for ( ii = 0; ii < 4; ii++ )
467 xx[cnt] = x[
items[ii].item];
468 yy[cnt] = y[
items[ii].item];
473 d1 = sqrt( ( xx[1] - xx[0] ) * ( xx[1] - xx[0] ) + ( yy[1] - yy[0] ) * ( yy[1] - yy[0] ) );
474 d2 = sqrt( ( xx[2] - xx[1] ) * ( xx[2] - xx[1] ) + ( yy[2] - yy[1] ) * ( yy[2] - yy[1] ) );
475 d3 = sqrt( ( xx[0] - xx[2] ) * ( xx[0] - xx[2] ) + ( yy[0] - yy[2] ) * ( yy[0] - yy[2] ) );
476 if ( d1 == 0. || d2 == 0. || d3 == 0. )
482 t = d1; d1 = d2; d2 = t;
487 t = d2; d2 = d3; d3 = t;
490 t = ( d1 + d2 ) / d3;
498 if ( excl_item == -1 )
503 for ( ii = 0; ii < 4; ii++ )
505 if ( ii != excl_item )
507 xx[cnt] = x[
items[ii].item];
508 yy[cnt] = y[
items[ii].item];
509 zz[cnt] = z[
items[ii].item];
514 A = yy[0] * ( zz[1] - zz[2] ) + yy[1] * ( zz[2] - zz[0] ) + yy[2] * ( zz[0] - zz[1] );
515 B = zz[0] * ( xx[1] - xx[2] ) + zz[1] * ( xx[2] - xx[0] ) + zz[2] * ( xx[0] - xx[1] );
516 C = xx[0] * ( yy[1] - yy[2] ) + xx[1] * ( yy[2] - yy[0] ) + xx[2] * ( yy[0] - yy[1] );
517 D = -A * xx[0] - B * yy[0] - C * zz[0];
520 zops->set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
541 for ( i = 0; i < nptsx; i++ )
543 for ( j = 0; j < nptsy; j++ )
545 dist2( xg[i], yg[j], x, y, npts );
546 zops->set( zgp, i, j, 0. );
548 for ( k = 0; k < 4; k++ )
550 if (
items[k].item != -1 )
553 zops->add( zgp, i, j, d * z[
items[k].item] );
558 zops->set( zgp, i, j,
NaN );
560 zops->div( zgp, i, j, nt );
587 if ( ( pin = (
point *) malloc( (
size_t) npts *
sizeof (
point ) ) ) == NULL )
589 plexit(
"grid_dtli: Insufficient memory" );
596 for ( i = 0; i < npts; i++ )
598 pt->x = (double) *xt++;
599 pt->y = (double) *yt++;
600 pt->z = (double) *zt++;
604 nptsg = nptsx * nptsy;
606 if ( ( pgrid = (
point *) malloc( (
size_t) nptsg *
sizeof (
point ) ) ) == NULL )
608 plexit(
"grid_dtli: Insufficient memory" );
613 for ( j = 0; j < nptsy; j++ )
616 for ( i = 0; i < nptsx; i++ )
618 pt->x = (double) *xt++;
619 pt->y = (double) *yt;
626 for ( i = 0; i < nptsx; i++ )
628 for ( j = 0; j < nptsy; j++ )
630 pt = &pgrid[j * nptsx + i];
631 zops->set( zgp, i, j, (
PLFLT)
pt->z );
660 plwarn(
"plgriddata(): GRID_NNI: wtmin must be specified with 'data' arg. Using -PLFLT_MAX" );
664 if ( ( pin = (
point *) malloc( (
size_t) npts *
sizeof (
point ) ) ) == NULL )
666 plexit(
"plgridata: Insufficient memory" );
673 for ( i = 0; i < npts; i++ )
675 pt->x = (double) *xt++;
676 pt->y = (double) *yt++;
677 pt->z = (double) *zt++;
681 nptsg = nptsx * nptsy;
683 if ( ( pgrid = (
point *) malloc( (
size_t) nptsg *
sizeof (
point ) ) ) == NULL )
685 plexit(
"plgridata: Insufficient memory" );
690 for ( j = 0; j < nptsy; j++ )
693 for ( i = 0; i < nptsx; i++ )
695 pt->x = (double) *xt++;
696 pt->y = (double) *yt;
703 for ( i = 0; i < nptsx; i++ )
705 for ( j = 0; j < nptsy; j++ )
707 pt = &pgrid[j * nptsx + i];
708 zops->set( zgp, i, j, (
PLFLT)
pt->z );
731 for ( i = 0; i < knn_order; i++ )
737 for ( i = 0; i < npts; i++ )
739 d = ( ( gx - x[i] ) * ( gx - x[i] ) + ( gy - y[i] ) * ( gy - y[i] ) );
747 items[max_slot].dist = d;
748 items[max_slot].item = i;
751 max_dist =
items[0].dist;
753 for ( j = 1; j < knn_order; j++ )
755 if (
items[j].dist > max_dist )
757 max_dist =
items[j].dist;
763 for ( j = 0; j < knn_order; j++ )
778 for ( i = 0; i < 4; i++ )
784 for ( i = 0; i < npts; i++ )
786 d = ( ( gx - x[i] ) * ( gx - x[i] ) + ( gy - y[i] ) * ( gy - y[i] ) );
792 quad = 2 * ( x[i] > gx ) + ( y[i] < gy );
798 if ( d <
items[quad].dist )
800 items[quad].dist = d;
801 items[quad].item = i;
805 for ( i = 0; i < 4; i++ )
806 if (
items[i].item != -1 )
817 boolT ismalloc = False;
820 vertexT *vertex, **vertexp;
821 facetT *neighbor, **neighborp;
822 int curlong, totlong;
823 FILE *outfile = NULL;
824 FILE *errfile = stderr;
829 PLFLT xt[3], yt[3], zt[3];
835 int numfacets, numsimplicial, numridges;
836 int totneighbors, numcoplanars, numtricoplanars;
838 plwarn(
"plgriddata: GRID_DTLI, If you have QHull knowledge, FIXME." );
842 strcpy( flags,
"qhull d Qbb Qt", 250 );
844 if ( ( points = (coordT *) malloc( npts * ( dim + 1 ) *
sizeof ( coordT ) ) ) == NULL )
846 plexit(
"grid_adtli: Insufficient memory" );
849 for ( i = 0; i < npts; i++ )
851 points[i * dim] = x[i];
852 points[i * dim + 1] = y[i];
856 exitcode = qh_new_qhull( dim, npts, points, ismalloc,
857 flags, outfile, errfile );
859 qh_init_A( stdin, stdout, stderr, 0, NULL );
860 exitcode = setjmp( qh errexit );
863 qh_initflags( flags );
864 qh PROJECTdelaunay = True;
865 qh_init_B( points, npts, dim, ismalloc );
873 printf(
"Triangles\n" );
875 if ( !facet->upperdelaunay )
877 FOREACHvertex_( facet->vertices )
878 printf( " %d", qh_pointid( vertex->
point ) );
885 printf(
"Neigbors\n" );
887 qh_findgood_all( qh facet_list );
888 qh_countfacets( qh facet_list, NULL, !qh_ALL, &numfacets, &numsimplicial,
889 &totneighbors, &numridges, &numcoplanars, &numtricoplanars );
892 if ( !facet->upperdelaunay )
894 FOREACHneighbor_( facet )
895 printf(
" %d", neighbor->visitid ? neighbor->visitid - 1 : -neighbor->id );
902 exitcode = setjmp( qh errexit );
905 qh NOerrexit = False;
906 for ( i = 0; i < nptsx; i++ )
907 for ( j = 0; j < nptsy; j++ )
912 qh_setdelaunay( 3, 1,
point );
918 facet = qh_findbestfacet(
point, qh_ALL, &bestdist, &isoutside );
922 facet = qh_findbest(
point, qh facet_list, qh_ALL,
925 &bestdist, &isoutside, &totpart );
929 vertex = qh_nearvertex( facet,
point, &bestdist );
945 facet = qh_findfacet_all(
point, &bestdist, &isoutside, &totpart );
947 if ( facet->upperdelaunay )
948 zops->set( zgp, i, j,
NaN );
951 FOREACHvertex_( facet->vertices )
953 k = qh_pointid( vertex->point );
962 A = yt[0] * ( zt[1] - zt[2] ) + yt[1] * ( zt[2] - zt[0] ) + yt[2] * ( zt[0] - zt[1] );
963 B = zt[0] * ( xt[1] - xt[2] ) + zt[1] * ( xt[2] - xt[0] ) + zt[2] * ( xt[0] - xt[1] );
964 C = xt[0] * ( yt[1] - yt[2] ) + xt[1] * ( yt[2] - yt[0] ) + xt[2] * ( yt[0] - yt[1] );
965 D = -A * xt[0] - B * yt[0] -
C * zt[0];
968 zops->set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
976 qh_freeqhull( !qh_ALL );
977 qh_memfreeshort( &curlong, &totlong );
978 if ( curlong || totlong )
980 "qhull: did not free %d bytes of long memory (%d pieces)\n",
void csa_addpoints(csa *a, int n, point points[])
void csa_calculatespline(csa *a)
void csa_approximate_points(csa *a, int n, point *points)
void lpi_interpolate_points(int nin, point pin[], int nout, point pout[])
NNDLLIMPEXP void nnpi_interpolate_points(int nin, point pin[], double wmin, int nout, point pout[])
void plwarn(PLCHAR_VECTOR errormsg)
void plexit(PLCHAR_VECTOR errormsg)
void plabort(PLCHAR_VECTOR errormsg)
void plfgriddata(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, PLINT npts, PLFLT_VECTOR xg, PLINT nptsx, PLFLT_VECTOR yg, PLINT nptsy, PLF2OPS zops, PLPointer zgp, PLINT type, PLFLT data)
static void dist2(PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts)
static void grid_nnli(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts, PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp, PLFLT threshold)
static PT items[KNN_MAX_ORDER]
void c_plgriddata(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, PLINT npts, PLFLT_VECTOR xg, PLINT nptsx, PLFLT_VECTOR yg, PLINT nptsy, PLFLT **zg, PLINT type, PLFLT data)
static void grid_nnidw(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts, PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp, int knn_order)
static void grid_nnaidw(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts, PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp)
static void dist1(PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts, int knn_order)
const PLFLT * PLFLT_VECTOR
plgriddata(x, y, z, xg, yg, type, data)\n\ \n\ \n\ This function is used in example 21.\n\ \n\ \n\ \n\ SYNOPSIS:\n\ \n\ plgriddata(x, y, z, npts, xg, nptsx, yg, nptsy, zg, type, data)\n\ \n\ ARGUMENTS:\n\ \n\ x(PLFLT_VECTOR, input) : The input x vector.\n\ \n\ y(PLFLT_VECTOR, input) : The input y vector.\n\ \n\ z(PLFLT_VECTOR, input) : The input z vector. Each triple x[i],\n\ y[i], z[i] represents one data sample coordinate.\n\ \n\ npts(PLINT, input) : The number of data samples in the x, y and z\n\ vectors.\n\ \n\ xg(PLFLT_VECTOR, input) : A vector that specifies the grid spacing\n\ in the x direction. Usually xg has nptsx equally spaced values\n\ from the minimum to the maximum values of the x input vector.\n\ \n\ nptsx(PLINT, input) : The number of points in the xg vector.\n\ \n\ yg(PLFLT_VECTOR, input) : A vector that specifies the grid spacing\n\ in the y direction. Similar to the xg parameter.\n\ \n\ nptsy(PLINT, input) : The number of points in the yg vector.\n\ \n\ zg(PLFLT_NC_MATRIX, output) : The matrix of interpolated results\n\ where data lies in the grid specified by xg and yg. Therefore the\n\ zg matrix must be dimensioned\n\ nptsx by\n\ nptsy.\n\ \n\ type(PLINT, input) : The type of grid interpolation algorithm to\n\ use, which can be:GRID_CSA:Bivariate Cubic Spline approximation\n\ GRID_DTLI:Delaunay Triangulation Linear Interpolation\n\ GRID_NNI:Natural Neighbors Interpolation\n\ GRID_NNIDW:Nearest Neighbors Inverse Distance Weighted\n\ GRID_NNLI:Nearest Neighbors Linear Interpolation\n\ GRID_NNAIDW: Nearest Neighbors Around Inverse Distance\n\ Weighted\n\ For details of the algorithms read the source file plgridd.c.\n\ \n\ data(PLFLT, input) : Some gridding algorithms require extra data,\n\ which can be specified through this argument. Currently, for\n\ algorithm:GRID_NNIDW, data specifies the number of neighbors to\n\ use, the lower the value, the noisier(more local) the\n\ approximation is.\n\ GRID_NNLI, data specifies what a thin triangle is, in the\n\ range[1. .. 2.]. High values enable the usage of very thin\n\ triangles for interpolation, possibly resulting in error in\n\ the approximation.\n\ GRID_NNI, only weights greater than data will be accepted. If\n\ 0, all weights will be accepted.\n\ " zg
static const char shade or gradient plots n n or n gradient plots(See pllegend for similar functionality for creating\n\ legends with discrete elements). The arguments of plcolorbar provide\n\ control over the location and size of the color bar as well as the\n\ location and characteristics of the elements(most of which are\n\ optional) within that color bar. The resulting color bar is clipped\n\ at the boundaries of the current subpage.(N.B. the adopted coordinate\n\ system used for some of the parameters is defined in the documentation\n\ of the position parameter.)\n\ \n\ Redacted form reads the desired grid location from the input vectors n xg[nptsx] and yg[nptsy]