// Grid data demo // // Copyright (C) 2009 Werner Smekal // // This file is part of PLplot. // // PLplot is free software; you can redistribute it and/or modify // it under the terms of the GNU Library General Public License as published // by the Free Software Foundation; either version 2 of the License, or // (at your option) any later version. // // PLplot is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Library General Public License for more details. // // You should have received a copy of the GNU Library General Public License // along with PLplot; if not, write to the Free Software // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA // // import std.math; import plplot; // Options data structure definition. PLINT pts = 500; PLINT xp = 25; PLINT yp = 20; PLINT nl = 16; int knn_order = 20; PLFLT threshold = 1.001; PLFLT wmin = -1e3; int randn = 0; int rosen = 0; PLFLT xm, xM, ym, yM; int main( char[][] args ) { string[] title = [ "Cubic Spline Approximation", "Delaunay Linear Interpolation", "Natural Neighbors Interpolation", "KNN Inv. Distance Weighted", "3NN Linear Interpolation", "4NN Around Inv. Dist. Weighted" ]; xm = ym = -0.2; xM = yM = 0.6; // plMergeOpts(options, "x21c options", NULL); plparseopts( args, PL_PARSE_FULL ); PLFLT[] opt = [ 0.0, 0.0, wmin, knn_order, threshold, 0.0 ]; // Initialize plplot plinit(); cmap1_init(); // Initialise random number generator plseed( 5489 ); PLFLT[] x, y, z; x.length = y.length = z.length = pts; create_data( x, y, z ); // the sampled data PLFLT zmin = z[0]; PLFLT zmax = z[0]; for ( int i = 1; i < pts; i++ ) { if ( z[i] > zmax ) zmax = z[i]; if ( z[i] < zmin ) zmin = z[i]; } PLFLT[] xg, yg; xg.length = xp; yg.length = yp; create_grid( xg, yg ); // grid the data at PLFLT[][] zg = new PLFLT[][xp]; for ( int i = 0; i < xp; i++ ) zg[i] = new PLFLT[yp]; PLFLT[] clev = new PLFLT[nl]; PLFLT[] xx = new PLFLT[1]; PLFLT[] yy = new PLFLT[1]; plcol0( 1 ); plenv( xm, xM, ym, yM, 2, 0 ); plcol0( 15 ); pllab( "X", "Y", "The original data sampling" ); for ( int i = 0; i < pts; i++ ) { plcol1( ( z[i] - zmin ) / ( zmax - zmin ) ); xx[0] = x[i]; yy[0] = y[i]; plstring( xx, yy, "#(727)" ); } pladv( 0 ); plssub( 3, 2 ); for ( int k = 0; k < 2; k++ ) { pladv( 0 ); for ( int alg = 1; alg < 7; alg++ ) { plgriddata( x, y, z, xg, yg, zg, alg, opt[alg - 1] ); // - CSA can generate NaNs (only interpolates?!). // - DTLI and NNI can generate NaNs for points outside the convex hull // of the data points. // - NNLI can generate NaNs if a sufficiently thick triangle is not found // // PLplot should be NaN/Inf aware, but changing it now is quite a job... // so, instead of not plotting the NaN regions, a weighted average over // the neighbors is done. // if ( alg == GRID_CSA || alg == GRID_DTLI || alg == GRID_NNLI || alg == GRID_NNI ) { PLFLT dist, d; for ( int i = 0; i < xp; i++ ) { for ( int j = 0; j < yp; j++ ) { if ( isNaN( zg[i][j] ) ) // average (IDW) over the 8 neighbors { zg[i][j] = 0.0; dist = 0.0; for ( int ii = i - 1; ii <= i + 1 && ii < xp; ii++ ) { for ( int jj = j - 1; jj <= j + 1 && jj < yp; jj++ ) { if ( ii >= 0 && jj >= 0 && !isNaN( zg[ii][jj] ) ) { d = ( abs( ii - i ) + abs( jj - j ) ) == 1 ? 1.0 : 1.4142; zg[i][j] += zg[ii][jj] / ( d * d ); dist += d; } } } if ( dist != 0.0 ) zg[i][j] /= dist; else zg[i][j] = zmin; } } } } PLFLT lzM, lzm; plMinMax2dGrid( zg, lzM, lzm ); lzm = fmin( lzm, zmin ); lzM = fmax( lzM, zmax ); // Increase limits slightly to prevent spurious contours // due to rounding errors lzm = lzm - 0.01; lzM = lzM + 0.01; plcol0( 1 ); pladv( alg ); if ( k == 0 ) { for ( int i = 0; i < nl; i++ ) clev[i] = lzm + ( lzM - lzm ) / ( nl - 1 ) * i; plenv0( xm, xM, ym, yM, 2, 0 ); plcol0( 15 ); pllab( "X", "Y", title[alg - 1] ); plshades( zg, null, xm, xM, ym, yM, clev, 1, 0, 1, 1 ); plcol0( 2 ); } else { for ( int i = 0; i < nl; i++ ) clev[i] = lzm + ( lzM - lzm ) / ( nl - 1 ) * i; plvpor( 0.0, 1.0, 0.0, 0.9 ); plwind( -1.1, 0.75, -0.65, 1.20 ); // // For the comparison to be fair, all plots should have the // same z values, but to get the max/min of the data generated // by all algorithms would imply two passes. Keep it simple. // // plw3d(1., 1., 1., xm, xM, ym, yM, zmin, zmax, 30, -60); // plw3d( 1., 1., 1., xm, xM, ym, yM, lzm, lzM, 30, -40 ); plbox3( "bntu", "X", 0., 0, "bntu", "Y", 0., 0, "bcdfntu", "Z", 0.5, 0 ); plcol0( 15 ); pllab( "", "", title[alg - 1] ); plot3dc( xg, yg, zg, DRAW_LINEXY | MAG_COLOR | BASE_CONT, clev ); } } } plend(); return 0; } void create_grid( PLFLT[] x, PLFLT[] y ) { int px = cast(int) x.length; int py = cast(int) y.length; for ( int i = 0; i < px; i++ ) x[i] = xm + ( xM - xm ) * i / ( px - 1.0 ); for ( int i = 0; i < py; i++ ) y[i] = ym + ( yM - ym ) * i / ( py - 1.0 ); } void create_data( PLFLT[] x, PLFLT[] y, PLFLT[] z ) { int pts = cast(int) x.length; assert( pts == y.length, "create_data(): Arrays must be of same length" ); assert( pts == z.length, "create_data(): Arrays must be of same length" ); PLFLT xt, yt, r; for ( int i = 0; i < pts; i++ ) { xt = ( xM - xm ) * plrandd(); yt = ( yM - ym ) * plrandd(); if ( !randn ) { x[i] = xt + xm; y[i] = yt + ym; } else // std=1, meaning that many points are outside the plot range { x[i] = sqrt( -2.0 * log( xt ) ) * cos( 2. * PI * yt ) + xm; y[i] = sqrt( -2.0 * log( xt ) ) * sin( 2. * PI * yt ) + ym; } if ( !rosen ) { r = sqrt( x[i] * x[i] + y[i] * y[i] ); z[i] = exp( -r * r ) * cos( 2.0 * PI * r ); } else z[i] = log( pow( 1. - x[i], 2.9 ) + 100.0 * pow( y[i] - pow( x[i], 2.0 ), 2.0 ) ); } } void cmap1_init() { PLFLT[] i = [ 0.0, 1.0 ]; // boundaries PLFLT[] h = [ 240.0, 0.0 ]; // blue -> green -> yellow -> red PLFLT[] l = [ 0.6, 0.6 ]; PLFLT[] s = [ 0.8, 0.8 ]; plscmap1n( 256 ); plscmap1l( 0, i, h, l, s ); }