1 // Grid data demo 2 // 3 // Copyright (C) 2004 Joao Cardoso 4 // 5 // This file is part of PLplot. 6 // 7 // PLplot is free software; you can redistribute it and/or modify 8 // it under the terms of the GNU Library General Public License as published 9 // by the Free Software Foundation; either version 2 of the License, or 10 // (at your option) any later version. 11 // 12 // PLplot is distributed in the hope that it will be useful, 13 // but WITHOUT ANY WARRANTY; without even the implied warranty of 14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 // GNU Library General Public License for more details. 16 // 17 // You should have received a copy of the GNU Library General Public License 18 // along with PLplot; if not, write to the Free Software 19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 20 // 21 // 22 23 #include "plcdemos.h" 24 25 // Options data structure definition. 26 27 static PLINT pts = 500; 28 static PLINT xp = 25; 29 static PLINT yp = 20; 30 static PLINT nl = 16; 31 static int knn_order = 20; 32 static PLFLT threshold = 1.001; 33 static PLFLT wmin = -1e3; 34 static int randn = 0; 35 static int rosen = 0; 36 37 static PLOptionTable options[] = { 38 { 39 "npts", 40 NULL, 41 NULL, 42 &pts, 43 PL_OPT_INT, 44 "-npts points", 45 "Specify number of random points to generate [500]" 46 }, 47 { 48 "randn", 49 NULL, 50 NULL, 51 &randn, 52 PL_OPT_BOOL, 53 "-randn", 54 "Normal instead of uniform sampling -- the effective \n\ 55 \t\t\t number of points will be smaller than the specified." 56 }, 57 { 58 "rosen", 59 NULL, 60 NULL, 61 &rosen, 62 PL_OPT_BOOL, 63 "-rosen", 64 "Generate points from the Rosenbrock function." 65 }, 66 { 67 "nx", 68 NULL, 69 NULL, 70 &xp, 71 PL_OPT_INT, 72 "-nx points", 73 "Specify grid x dimension [25]" 74 }, 75 { 76 "ny", 77 NULL, mypltr(PLFLT x,PLFLT y,PLFLT * tx,PLFLT * ty,PLPointer pltr_data)78 NULL, 79 &yp, 80 PL_OPT_INT, 81 "-ny points", 82 "Specify grid y dimension [20]" 83 }, 84 { 85 "nlevel", 86 NULL, 87 NULL, 88 &nl, 89 PL_OPT_INT, 90 "-nlevel ", 91 "Specify number of contour levels [15]" 92 }, 93 { 94 "knn_order", 95 NULL, 96 NULL, 97 &knn_order, 98 PL_OPT_INT, 99 "-knn_order order", 100 "Specify the number of neighbors [20]" 101 }, 102 { 103 "threshold", 104 NULL, 105 NULL, 106 &threshold, 107 PL_OPT_FLOAT, 108 "-threshold float", 109 "Specify what a thin triangle is [1. < [1.001] < 2.]" 110 }, 111 112 { 113 NULL, // option 114 NULL, // handler 115 NULL, // client data 116 NULL, // address of variable to set 117 0, // mode flag 118 NULL, // short syntax 119 NULL 120 } // long syntax 121 }; 122 123 void create_data( PLFLT **xi, PLFLT **yi, PLFLT **zi, int npts ); 124 void free_data( PLFLT *x, PLFLT *y, PLFLT *z ); 125 void create_grid( PLFLT **xi, int px, PLFLT **yi, int py ); 126 void free_grid( PLFLT *x, PLFLT *y ); 127 128 static void 129 cmap1_init( void ) 130 { 131 PLFLT i[2], h[2], l[2], s[2]; 132 133 i[0] = 0.0; // left boundary 134 i[1] = 1.0; // right boundary 135 136 h[0] = 240; // blue -> green -> yellow -> 137 h[1] = 0; // -> red 138 139 l[0] = 0.6; 140 l[1] = 0.6; 141 142 s[0] = 0.8; 143 s[1] = 0.8; 144 145 plscmap1n( 256 ); 146 c_plscmap1l( 0, 2, i, h, l, s, NULL ); 147 } 148 149 PLFLT xm, xM, ym, yM; 150 151 int 152 main( int argc, char *argv[] ) 153 { 154 PLFLT *x, *y, *z, *clev; 155 PLFLT *xg, *yg, **zg; 156 PLFLT zmin, zmax, lzm, lzM; 157 int i, j, k; 158 PLINT alg; 159 PLCHAR_VECTOR title[] = { "Cubic Spline Approximation", 160 "Delaunay Linear Interpolation", 161 "Natural Neighbors Interpolation", 162 "KNN Inv. Distance Weighted", 163 "3NN Linear Interpolation", 164 "4NN Around Inv. Dist. Weighted" }; 165 166 PLFLT opt[] = { 0., 0., 0., 0., 0., 0. }; 167 168 xm = ym = -0.2; 169 xM = yM = 0.6; 170 171 plMergeOpts( options, "x21c options", NULL ); 172 plparseopts( &argc, argv, PL_PARSE_FULL ); 173 174 opt[2] = wmin; 175 opt[3] = (PLFLT) knn_order; 176 opt[4] = threshold; 177 178 // Initialize plplot 179 180 plinit(); 181 182 // Use a colour map with no black band in the middle. 183 cmap1_init(); 184 // Initialise random number generator 185 plseed( 5489 ); 186 187 create_data( &x, &y, &z, pts ); // the sampled data 188 zmin = z[0]; 189 zmax = z[0]; 190 for ( i = 1; i < pts; i++ ) 191 { 192 if ( z[i] > zmax ) 193 zmax = z[i]; 194 if ( z[i] < zmin ) 195 zmin = z[i]; 196 } 197 198 create_grid( &xg, xp, &yg, yp ); // grid the data at 199 plAlloc2dGrid( &zg, xp, yp ); // the output grided data 200 clev = (PLFLT *) malloc( (size_t) nl * sizeof ( PLFLT ) ); 201 202 plcol0( 1 ); 203 plenv( xm, xM, ym, yM, 2, 0 ); 204 plcol0( 15 ); 205 pllab( "X", "Y", "The original data sampling" ); 206 for ( i = 0; i < pts; i++ ) 207 { 208 plcol1( ( z[i] - zmin ) / ( zmax - zmin ) ); 209 // The following plstring call should be the the equivalent of 210 // plpoin( 1, &x[i], &y[i], 5 ); Use plstring because it is 211 // not deprecated like plpoin and has much more powerful 212 // capabilities. N.B. symbol 141 works for Hershey devices 213 // (e.g., -dev xwin) only if plfontld( 0 ) has been called 214 // while symbol 727 works only if plfontld( 1 ) has been 215 // called. The latter is the default which is why we use 727 216 // here to represent a centred X (multiplication) symbol. 217 // This dependence on plfontld is one of the limitations of 218 // the Hershey escapes for PLplot, but the upside is you get 219 // reasonable results for both Hershey and Unicode devices. 220 plstring( 1, &x[i], &y[i], "#(727)" ); 221 } 222 pladv( 0 ); 223 224 plssub( 3, 2 ); 225 226 for ( k = 0; k < 2; k++ ) 227 { 228 pladv( 0 ); 229 for ( alg = 1; alg < 7; alg++ ) 230 { 231 plgriddata( x, y, z, pts, xg, xp, yg, yp, zg, alg, opt[alg - 1] ); 232 233 // - CSA can generate NaNs (only interpolates?!). 234 // - DTLI and NNI can generate NaNs for points outside the convex hull 235 // of the data points. 236 // - NNLI can generate NaNs if a sufficiently thick triangle is not found 237 // 238 // PLplot should be NaN/Inf aware, but changing it now is quite a job... 239 // so, instead of not plotting the NaN regions, a weighted average over 240 // the neighbors is done. 241 // 242 243 if ( alg == GRID_CSA || alg == GRID_DTLI || alg == GRID_NNLI || alg == GRID_NNI ) 244 { 245 int ii, jj; 246 PLFLT dist, d; 247 248 for ( i = 0; i < xp; i++ ) 249 { 250 for ( j = 0; j < yp; j++ ) 251 { 252 if ( isnan( zg[i][j] ) ) // average (IDW) over the 8 neighbors 253 254 { 255 zg[i][j] = 0.; dist = 0.; 256 257 for ( ii = i - 1; ii <= i + 1 && ii < xp; ii++ ) 258 { 259 for ( jj = j - 1; jj <= j + 1 && jj < yp; jj++ ) 260 { 261 if ( ii >= 0 && jj >= 0 && !isnan( zg[ii][jj] ) ) 262 { 263 d = ( abs( ii - i ) + abs( jj - j ) ) == 1 ? 1. : 1.4142; 264 zg[i][j] += zg[ii][jj] / ( d * d ); 265 dist += d; 266 } 267 } 268 } 269 if ( dist != 0. ) 270 zg[i][j] /= dist; 271 else 272 zg[i][j] = zmin; 273 } 274 } 275 } 276 } 277 278 plMinMax2dGrid( (PLFLT_MATRIX) zg, xp, yp, &lzM, &lzm ); 279 280 lzm = MIN( lzm, zmin ); 281 lzM = MAX( lzM, zmax ); 282 283 // Increase limits slightly to prevent spurious contours 284 // due to rounding errors 285 lzm = lzm - 0.01; 286 lzM = lzM + 0.01; 287 288 plcol0( 1 ); 289 290 pladv( alg ); 291 292 if ( k == 0 ) 293 { 294 for ( i = 0; i < nl; i++ ) 295 clev[i] = lzm + ( lzM - lzm ) / ( nl - 1 ) * i; 296 297 plenv0( xm, xM, ym, yM, 2, 0 ); 298 plcol0( 15 ); 299 pllab( "X", "Y", title[alg - 1] ); 300 plshades( (PLFLT_MATRIX) zg, xp, yp, NULL, xm, xM, ym, yM, 301 clev, nl, 1., 0, 1., plfill, 1, NULL, NULL ); 302 plcol0( 2 ); 303 } 304 else 305 { 306 for ( i = 0; i < nl; i++ ) 307 clev[i] = lzm + ( lzM - lzm ) / ( nl - 1 ) * i; 308 309 plvpor( 0.0, 1.0, 0.0, 0.9 ); 310 plwind( -1.1, 0.75, -0.65, 1.20 ); 311 // 312 // For the comparison to be fair, all plots should have the read_img(PLCHAR_VECTOR fname,PLFLT *** img_f,int * width,int * height,int * num_col)313 // same z values, but to get the max/min of the data generated 314 // by all algorithms would imply two passes. Keep it simple. 315 // 316 // plw3d(1., 1., 1., xm, xM, ym, yM, zmin, zmax, 30, -60); 317 // 318 319 plw3d( 1., 1., 1., xm, xM, ym, yM, lzm, lzM, 30, -40 ); 320 plbox3( "bntu", "X", 0., 0, 321 "bntu", "Y", 0., 0, 322 "bcdfntu", "Z", 0.5, 0 ); 323 plcol0( 15 ); 324 pllab( "", "", title[alg - 1] ); 325 plot3dc( xg, yg, (PLFLT_MATRIX) zg, xp, yp, DRAW_LINEXY | MAG_COLOR | BASE_CONT, clev, nl ); 326 } 327 } 328 } 329 330 plend(); 331 332 free_data( x, y, z ); 333 free_grid( xg, yg ); 334 free( (void *) clev ); 335 plFree2dGrid( zg, xp, yp ); 336 exit( 0 ); 337 } 338 339 340 void 341 create_grid( PLFLT **xi, int px, PLFLT **yi, int py ) 342 { 343 PLFLT *x, *y; 344 int i; 345 346 x = *xi = (PLFLT *) malloc( (size_t) px * sizeof ( PLFLT ) ); 347 y = *yi = (PLFLT *) malloc( (size_t) py * sizeof ( PLFLT ) ); 348 349 for ( i = 0; i < px; i++ ) 350 *x++ = xm + ( xM - xm ) * i / ( px - 1. ); 351 352 for ( i = 0; i < py; i++ ) 353 *y++ = ym + ( yM - ym ) * i / ( py - 1. ); 354 } 355 356 void 357 free_grid( PLFLT *xi, PLFLT *yi ) 358 { 359 free( (void *) xi ); 360 free( (void *) yi ); 361 } 362 363 void 364 create_data( PLFLT **xi, PLFLT **yi, PLFLT **zi, int npts ) 365 { 366 int i; 367 PLFLT *x, *y, *z, r; 368 PLFLT xt, yt; 369 370 *xi = x = (PLFLT *) malloc( (size_t) npts * sizeof ( PLFLT ) ); 371 *yi = y = (PLFLT *) malloc( (size_t) npts * sizeof ( PLFLT ) ); 372 *zi = z = (PLFLT *) malloc( (size_t) npts * sizeof ( PLFLT ) ); 373 374 for ( i = 0; i < npts; i++ ) 375 { 376 xt = ( xM - xm ) * plrandd(); 377 yt = ( yM - ym ) * plrandd(); 378 if ( !randn ) 379 { 380 *x = xt + xm; 381 *y = yt + ym; save_plot(char * fname)382 } 383 else // std=1, meaning that many points are outside the plot range 384 { 385 *x = sqrt( -2. * log( xt ) ) * cos( 2. * M_PI * yt ) + xm; 386 *y = sqrt( -2. * log( xt ) ) * sin( 2. * M_PI * yt ) + ym; 387 } 388 if ( !rosen ) 389 { 390 r = sqrt( ( *x ) * ( *x ) + ( *y ) * ( *y ) ); 391 *z = exp( -r * r ) * cos( 2.0 * M_PI * r ); 392 } 393 else 394 { 395 *z = log( pow( 1. - *x, 2. ) + 100. * pow( *y - pow( *x, 2. ), 2. ) ); 396 } 397 x++; y++; z++; 398 } 399 } get_clip(PLFLT * xi,PLFLT * xe,PLFLT * yi,PLFLT * ye)400 401 void 402 free_data( PLFLT *x, PLFLT *y, PLFLT *z ) 403 { 404 free( (void *) x ); 405 free( (void *) y ); 406 free( (void *) z ); 407 } 408