29 #include "../lib/csa/csa.h"
31 #include "../lib/csa/nan.h"
34 #include "../lib/nn/nn.h"
35 #include <qhull/qhull_a.h>
78 #define KNN_MAX_ORDER 100
118 plfgriddata( x, y, z, npts, xg, nptsx, yg, nptsy,
plf2ops_c(), (
PLPointer) zg, type, data );
128 if ( npts < 1 || nptsx < 1 || nptsy < 1 )
130 plabort(
"plgriddata: Bad array dimensions" );
136 for ( i = 0; i < nptsx - 1; i++ )
138 if ( xg[i] >= xg[i + 1] )
140 plabort(
"plgriddata: xg array must be strictly increasing" );
144 for ( i = 0; i < nptsy - 1; i++ )
146 if ( yg[i] >= yg[i + 1] )
148 plabort(
"plgriddata: yg array must be strictly increasing" );
154 for ( i = 0; i < nptsx; i++ )
155 for ( j = 0; j < nptsy; j++ )
156 zops->
set( zgp, i, j, 0.0 );
163 grid_csa( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
165 plwarn(
"plgriddata(): PLplot was configured to not use GRID_CSA.\n Reverting to GRID_NNAIDW." );
166 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
171 grid_nnidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, (
int) data );
175 grid_nnli( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, data );
179 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
184 grid_dtli( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
186 plwarn(
"plgriddata(): you must have the Qhull library installed to use GRID_DTLI.\n Reverting to GRID_NNAIDW." );
187 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
193 grid_nni( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, data );
195 plwarn(
"plgriddata(): you must have the Qhull library installed to use GRID_NNI.\n Reverting to GRID_NNAIDW." );
196 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
201 plabort(
"plgriddata: unknown algorithm type" );
217 const PLFLT *xt, *yt, *zt;
222 if ( ( pin = (
point *) malloc( (
size_t) npts *
sizeof (
point ) ) ) == NULL )
224 plexit(
"grid_csa: Insufficient memory" );
231 for ( i = 0; i < npts; i++ )
233 pt->
x = (double) *xt++;
234 pt->
y = (double) *yt++;
235 pt->
z = (double) *zt++;
239 nptsg = nptsx * nptsy;
240 if ( ( pgrid = (
point *) malloc( (
size_t) nptsg *
sizeof (
point ) ) ) == NULL )
242 plexit(
"grid_csa: Insufficient memory" );
247 for ( j = 0; j < nptsy; j++ )
250 for ( i = 0; i < nptsx; i++ )
252 pt->
x = (double) *xt++;
253 pt->
y = (double) *yt;
264 for ( i = 0; i < nptsx; i++ )
266 for ( j = 0; j < nptsy; j++ )
268 pt = &pgrid[j * nptsx + i];
289 const PLFLT *xg,
int nptsx,
const PLFLT *yg,
int nptsy,
297 plabort(
"plgriddata(): GRID_NNIDW: knn_order too big" );
301 if ( knn_order == 0 )
303 plwarn(
"plgriddata(): GRID_NNIDW: knn_order must be specified with 'data' arg. Using 15" );
307 for ( i = 0; i < nptsx; i++ )
309 for ( j = 0; j < nptsy; j++ )
311 dist1( xg[i], yg[j], x, y, npts, knn_order );
313 #ifdef GMS // alternative weight coeficients. I Don't like the results
316 for ( k = 1; k < knn_order; k++ )
317 if ( items[k].dist > md )
320 zops->
set( zgp, i, j, 0.0 );
323 for ( k = 0; k < knn_order; k++ )
325 if ( items[k].item == -1 )
328 wi = ( md - items[k].
dist ) / ( md * items[k].dist );
331 wi = 1. / ( items[k].
dist * items[k].
dist );
333 zops->
add( zgp, i, j, wi * z[items[k].item] );
337 zops->
div( zgp, i, j, nt );
339 zops->
set( zgp, i, j,
NaN );
352 const PLFLT *xg,
int nptsx,
const PLFLT *yg,
int nptsy,
355 PLFLT xx[4], yy[4], zz[4], t, A, B, C, D, d1, d2, d3, max_thick;
356 int i, j, ii, excl, cnt, excl_item;
358 if ( threshold == 0. )
360 plwarn(
"plgriddata(): GRID_NNLI: threshold must be specified with 'data' arg. Using 1.001" );
363 else if ( threshold > 2. || threshold < 1. )
365 plabort(
"plgriddata(): GRID_NNLI: 1. < threshold < 2." );
369 for ( i = 0; i < nptsx; i++ )
371 for ( j = 0; j < nptsy; j++ )
373 dist1( xg[i], yg[j], x, y, npts, 3 );
376 for ( ii = 0; ii < 3; ii++ )
378 xx[ii] = x[items[ii].
item];
379 yy[ii] = y[items[ii].item];
380 zz[ii] = z[items[ii].item];
383 d1 = sqrt( ( xx[1] - xx[0] ) * ( xx[1] - xx[0] ) + ( yy[1] - yy[0] ) * ( yy[1] - yy[0] ) );
384 d2 = sqrt( ( xx[2] - xx[1] ) * ( xx[2] - xx[1] ) + ( yy[2] - yy[1] ) * ( yy[2] - yy[1] ) );
385 d3 = sqrt( ( xx[0] - xx[2] ) * ( xx[0] - xx[2] ) + ( yy[0] - yy[2] ) * ( yy[0] - yy[2] ) );
387 if ( d1 == 0. || d2 == 0. || d3 == 0. )
389 zops->
set( zgp, i, j,
NaN );
396 t = d1; d1 = d2; d2 = t;
402 t = d2; d2 = d3; d3 = t;
405 if ( ( d1 + d2 ) / d3 < threshold )
407 zops->
set( zgp, i, j,
NaN );
412 A = yy[0] * ( zz[1] - zz[2] ) + yy[1] * ( zz[2] - zz[0] ) + yy[2] * ( zz[0] - zz[1] );
413 B = zz[0] * ( xx[1] - xx[2] ) + zz[1] * ( xx[2] - xx[0] ) + zz[2] * ( xx[0] - xx[1] );
414 C = xx[0] * ( yy[1] - yy[2] ) + xx[1] * ( yy[2] - yy[0] ) + xx[2] * ( yy[0] - yy[1] );
415 D = -A * xx[0] - B * yy[0] - C * zz[0];
418 zops->
set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
433 for ( i = 0; i < nptsx; i++ )
435 for ( j = 0; j < nptsy; j++ )
437 if ( zops->
is_nan( zgp, i, j ) )
439 dist1( xg[i], yg[j], x, y, npts, 4 );
453 max_thick = 0.; excl_item = -1;
454 for ( excl = 0; excl < 4; excl++ )
458 for ( ii = 0; ii < 4; ii++ )
462 xx[cnt] = x[items[ii].
item];
463 yy[cnt] = y[items[ii].item];
468 d1 = sqrt( ( xx[1] - xx[0] ) * ( xx[1] - xx[0] ) + ( yy[1] - yy[0] ) * ( yy[1] - yy[0] ) );
469 d2 = sqrt( ( xx[2] - xx[1] ) * ( xx[2] - xx[1] ) + ( yy[2] - yy[1] ) * ( yy[2] - yy[1] ) );
470 d3 = sqrt( ( xx[0] - xx[2] ) * ( xx[0] - xx[2] ) + ( yy[0] - yy[2] ) * ( yy[0] - yy[2] ) );
471 if ( d1 == 0. || d2 == 0. || d3 == 0. )
477 t = d1; d1 = d2; d2 = t;
482 t = d2; d2 = d3; d3 = t;
485 t = ( d1 + d2 ) / d3;
493 if ( excl_item == -1 )
498 for ( ii = 0; ii < 4; ii++ )
500 if ( ii != excl_item )
502 xx[cnt] = x[items[ii].
item];
503 yy[cnt] = y[items[ii].item];
504 zz[cnt] = z[items[ii].item];
509 A = yy[0] * ( zz[1] - zz[2] ) + yy[1] * ( zz[2] - zz[0] ) + yy[2] * ( zz[0] - zz[1] );
510 B = zz[0] * ( xx[1] - xx[2] ) + zz[1] * ( xx[2] - xx[0] ) + zz[2] * ( xx[0] - xx[1] );
511 C = xx[0] * ( yy[1] - yy[2] ) + xx[1] * ( yy[2] - yy[0] ) + xx[2] * ( yy[0] - yy[1] );
512 D = -A * xx[0] - B * yy[0] - C * zz[0];
515 zops->
set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
536 for ( i = 0; i < nptsx; i++ )
538 for ( j = 0; j < nptsy; j++ )
540 dist2( xg[i], yg[j], x, y, npts );
541 zops->
set( zgp, i, j, 0. );
543 for ( k = 0; k < 4; k++ )
545 if ( items[k].item != -1 )
547 d = 1. / ( items[k].
dist * items[k].
dist );
548 zops->
add( zgp, i, j, d * z[items[k].item] );
553 zops->
set( zgp, i, j,
NaN );
555 zops->
div( zgp, i, j, nt );
576 point *pin, *pgrid, *pt;
577 const PLFLT *xt, *yt, *zt;
580 if (
sizeof ( realT ) !=
sizeof (
double ) )
582 plabort(
"plgridata: QHull was compiled for floats instead of doubles" );
586 if ( ( pin = (
point *) malloc( (
size_t) npts *
sizeof (
point ) ) ) == NULL )
588 plexit(
"grid_dtli: Insufficient memory" );
595 for ( i = 0; i < npts; i++ )
597 pt->
x = (double) *xt++;
598 pt->
y = (double) *yt++;
599 pt->
z = (double) *zt++;
603 nptsg = nptsx * nptsy;
605 if ( ( pgrid = (
point *) malloc( (
size_t) nptsg *
sizeof (
point ) ) ) == NULL )
607 plexit(
"grid_dtli: Insufficient memory" );
612 for ( j = 0; j < nptsy; j++ )
615 for ( i = 0; i < nptsx; i++ )
617 pt->
x = (double) *xt++;
618 pt->
y = (double) *yt;
625 for ( i = 0; i < nptsx; i++ )
627 for ( j = 0; j < nptsy; j++ )
629 pt = &pgrid[j * nptsx + i];
650 const PLFLT *xt, *yt, *zt;
651 point *pin, *pgrid, *pt;
655 if (
sizeof ( realT ) !=
sizeof (
double ) )
657 plabort(
"plgridata: QHull was compiled for floats instead of doubles" );
663 plwarn(
"plgriddata(): GRID_NNI: wtmin must be specified with 'data' arg. Using -PLFLT_MAX" );
667 if ( ( pin = (
point *) malloc( (
size_t) npts *
sizeof (
point ) ) ) == NULL )
669 plexit(
"plgridata: Insufficient memory" );
676 for ( i = 0; i < npts; i++ )
678 pt->
x = (double) *xt++;
679 pt->
y = (double) *yt++;
680 pt->
z = (double) *zt++;
684 nptsg = nptsx * nptsy;
686 if ( ( pgrid = (
point *) malloc( (
size_t) nptsg *
sizeof (
point ) ) ) == NULL )
688 plexit(
"plgridata: Insufficient memory" );
693 for ( j = 0; j < nptsy; j++ )
696 for ( i = 0; i < nptsx; i++ )
698 pt->
x = (double) *xt++;
699 pt->
y = (double) *yt;
706 for ( i = 0; i < nptsx; i++ )
708 for ( j = 0; j < nptsy; j++ )
710 pt = &pgrid[j * nptsx + i];
718 #endif // PL_HAVE_QHULL
734 for ( i = 0; i < knn_order; i++ )
740 for ( i = 0; i < npts; i++ )
742 d = ( ( gx - x[i] ) * ( gx - x[i] ) + ( gy - y[i] ) * ( gy - y[i] ) );
750 items[max_slot].
dist = d;
751 items[max_slot].
item = i;
754 max_dist = items[0].
dist;
756 for ( j = 1; j < knn_order; j++ )
758 if ( items[j].dist > max_dist )
760 max_dist = items[j].
dist;
766 for ( j = 0; j < knn_order; j++ )
767 items[j].dist = sqrt( items[j].dist );
781 for ( i = 0; i < 4; i++ )
787 for ( i = 0; i < npts; i++ )
789 d = ( ( gx - x[i] ) * ( gx - x[i] ) + ( gy - y[i] ) * ( gy - y[i] ) );
795 quad = 2 * ( x[i] > gx ) + ( y[i] < gy );
801 if ( d < items[quad].dist )
803 items[quad].
dist = d;
804 items[quad].
item = i;
808 for ( i = 0; i < 4; i++ )
809 if ( items[i].item != -1 )
810 items[i].
dist = sqrt( items[i].dist );
814 #ifdef NONN // another DTLI, based only on QHULL, not nn
820 boolT ismalloc = False;
823 vertexT *vertex, **vertexp;
824 facetT *neighbor, **neighborp;
825 int curlong, totlong;
826 FILE *outfile = NULL;
827 FILE *errfile = stderr;
832 PLFLT xt[3], yt[3], zt[3];
838 int numfacets, numsimplicial, numridges;
839 int totneighbors, numcoplanars, numtricoplanars;
841 plwarn(
"plgriddata: GRID_DTLI, If you have QHull knowledge, FIXME." );
845 strcpy( flags,
"qhull d Qbb Qt", 250 );
847 if ( ( points = (coordT *) malloc( npts * ( dim + 1 ) *
sizeof ( coordT ) ) ) == NULL )
849 plexit(
"grid_adtli: Insufficient memory" );
852 for ( i = 0; i < npts; i++ )
854 points[i * dim] = x[i];
855 points[i * dim + 1] = y[i];
859 exitcode = qh_new_qhull( dim, npts, points, ismalloc,
860 flags, outfile, errfile );
862 qh_init_A( stdin, stdout, stderr, 0, NULL );
863 exitcode = setjmp( qh errexit );
866 qh_initflags( flags );
867 qh PROJECTdelaunay = True;
868 qh_init_B( points, npts, dim, ismalloc );
875 #if 0 // print the triangles vertices
876 printf(
"Triangles\n" );
878 if ( !facet->upperdelaunay )
880 FOREACHvertex_( facet->vertices )
881 printf( " %d", qh_pointid( vertex->point ) );
887 #if 0 // print each triangle neighbors
888 printf(
"Neigbors\n" );
890 qh_findgood_all( qh facet_list );
891 qh_countfacets( qh facet_list, NULL, !qh_ALL, &numfacets, &numsimplicial,
892 &totneighbors, &numridges, &numcoplanars, &numtricoplanars );
895 if ( !facet->upperdelaunay )
897 FOREACHneighbor_( facet )
898 printf( " %d", neighbor->visitid ? neighbor->visitid - 1 : -neighbor->
id );
905 exitcode = setjmp( qh errexit );
908 qh NOerrexit = False;
909 for ( i = 0; i < nptsx; i++ )
910 for ( j = 0; j < nptsy; j++ )
915 qh_setdelaunay( 3, 1, point );
921 facet = qh_findbestfacet( point, qh_ALL, &bestdist, &isoutside );
925 facet = qh_findbest( point, qh facet_list, qh_ALL,
928 &bestdist, &isoutside, &totpart );
932 vertex = qh_nearvertex( facet, point, &bestdist );
948 facet = qh_findfacet_all( point, &bestdist, &isoutside, &totpart );
950 if ( facet->upperdelaunay )
951 zops->
set( zgp, i, j,
NaN );
954 FOREACHvertex_( facet->vertices )
956 k = qh_pointid( vertex->point );
965 A = yt[0] * ( zt[1] - zt[2] ) + yt[1] * ( zt[2] - zt[0] ) + yt[2] * ( zt[0] - zt[1] );
966 B = zt[0] * ( xt[1] - xt[2] ) + zt[1] * ( xt[2] - xt[0] ) + zt[2] * ( xt[0] - xt[1] );
967 C = xt[0] * ( yt[1] - yt[2] ) + xt[1] * ( yt[2] - yt[0] ) + xt[2] * ( yt[0] - yt[1] );
968 D = -A * xt[0] - B * yt[0] - C * zt[0];
971 zops->
set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
979 qh_freeqhull( !qh_ALL );
980 qh_memfreeshort( &curlong, &totlong );
981 if ( curlong || totlong )
983 "qhull: did not free %d bytes of long memory (%d pieces)\n",
static void dist2(PLFLT gx, PLFLT gy, const PLFLT *x, const PLFLT *y, int npts)
PLFLT(* div)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
void plfgriddata(const PLFLT *x, const PLFLT *y, const PLFLT *z, PLINT npts, const PLFLT *xg, PLINT nptsx, const PLFLT *yg, PLINT nptsy, PLF2OPS zops, PLPointer zgp, PLINT type, PLFLT data)
void PLFLT PLINT PLINT PLFLT x
static PT items[KNN_MAX_ORDER]
void PLFLT PLINT PLINT PLFLT PLFLT PLFLT PLFLT PLINT PLINT PLINT PLFLT PLFLT PLINT PLFLT PLINT const PLINT const char *const PLINT const char *const const PLFLT const PLINT const PLINT const PLFLT * a
void csa_calculatespline(csa *a)
static void grid_nnaidw(const PLFLT *x, const PLFLT *y, const PLFLT *z, int npts, const PLFLT *xg, int nptsx, const PLFLT *yg, int nptsy, PLF2OPS zops, PLPointer zgp)
void csa_addpoints(csa *a, int n, point points[])
PLFLT(* set)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
void lpi_interpolate_points(int nin, point pin[], int nout, point pout[])
void c_plgriddata(const PLFLT *x, const PLFLT *y, const PLFLT *z, PLINT npts, const PLFLT *xg, PLINT nptsx, const PLFLT *yg, PLINT nptsy, PLFLT **zg, PLINT type, PLFLT data)
static void grid_nnli(const PLFLT *x, const PLFLT *y, const PLFLT *z, int npts, const PLFLT *xg, int nptsx, const PLFLT *yg, int nptsy, PLF2OPS zops, PLPointer zgp, PLFLT threshold)
void PLFLT PLINT PLINT PLFLT PLFLT y
static void dist1(PLFLT gx, PLFLT gy, const PLFLT *x, const PLFLT *y, int npts, int knn_order)
static void grid_nnidw(const PLFLT *x, const PLFLT *y, const PLFLT *z, int npts, const PLFLT *xg, int nptsx, const PLFLT *yg, int nptsy, PLF2OPS zops, PLPointer zgp, int knn_order)
void plabort(const char *errormsg)
void plwarn(const char *errormsg)
PLFLT(* add)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
NNDLLIMPEXP void nnpi_interpolate_points(int nin, point pin[], double wmin, int nout, point pout[])
void csa_approximate_points(csa *a, int n, point *points)
PLINT(* is_nan)(PLPointer p, PLINT ix, PLINT iy)
void plexit(const char *errormsg)