27#include "../lib/csa/csa.h"
29#include "../lib/csa/nan.h"
32#include "../lib/nn/nn.h"
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];
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 )
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++ )
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++ )
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 )
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];
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];
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] ) );
753 for ( j = 1; j < knn_order; j++ )
755 if (
items[j].dist > max_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 )
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
PLINT(* is_nan)(PLPointer p, PLINT ix, PLINT iy)
PLFLT(* add)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
PLFLT(* div)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
PLFLT(* set)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)