27 #include "../lib/csa/csa.h" 29 #include "../lib/csa/nan.h" 32 #include "../lib/nn/nn.h" 34 #ifdef PLPLOT_NONN // another DTLI, based only on QHULL, not nn 35 #ifdef HAS_LIBQHULL_INCLUDE 36 #include <libqhull/qhull_a.h> 37 #else //#ifdef HAS_LIBQHULL_INCLUDE 38 #include <qhull/qhull_a.h> 39 #endif //#ifdef HAS_LIBQHULL_INCLUDE 40 #endif //#ifdef PLPLOT_NONN 41 #endif //#ifdef PL_HAVE_QHUL 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 );
318 #ifdef GMS // alternative weight coeficients. I Don't like the results 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 );
336 wi = 1. / ( items[k].
dist * 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 )
552 d = 1. / ( items[k].
dist * items[k].
dist );
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];
715 #endif // PL_HAVE_QHULL 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++ )
764 items[j].
dist = sqrt( items[j].
dist );
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 )
807 items[i].
dist = sqrt( items[i].
dist );
811 #ifdef PLPLOT_NONN // another DTLI, based only on QHULL, not nn 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 );
872 #if 0 // print the triangles vertices 873 printf(
"Triangles\n" );
875 if ( !facet->upperdelaunay )
877 FOREACHvertex_( facet->vertices )
878 printf(
" %d", qh_pointid( vertex->point ) );
884 #if 0 // print each triangle neighbors 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",
983 #endif // PLPLOT_NONN
void plexit(PLCHAR_VECTOR errormsg)
void plabort(PLCHAR_VECTOR errormsg)
static PT items[KNN_MAX_ORDER]
PLFLT(* set)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
void csa_calculatespline(csa *a)
void csa_addpoints(csa *a, int n, point points[])
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 void dist2(PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts)
void lpi_interpolate_points(int nin, point pin[], int nout, point pout[])
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)
PLFLT(* add)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
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)
PLINT(* is_nan)(PLPointer p, PLINT ix, PLINT iy)
PLFLT(* div)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
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)
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)
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)
void plwarn(PLCHAR_VECTOR errormsg)
const PLFLT * PLFLT_VECTOR