libplplot-dev 5.10.0+dfsg-1 / usr / share / doc / libplplot12 / examples / d / x21d.d

// $Id: x21d.d 12393 2013-06-27 19:49:07Z andrewross $
//      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 );
}