#! /usr/bin/env perl
#
# Demo x21 for the PLplot PDL binding
#
# Grid data demo
#
# Copyright (C) 2004 Rafael Laboissiere
#
# 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
# SYNC: x21c.c 1.12
use PDL;
use PDL::Graphics::PLplot;
use Math::Trig qw [pi];
use Getopt::Long qw [:config pass_through];
use POSIX qw [clock];
sub cmap1_init {
my $i = pdl ([0.0, # left boundary
1.0]); # right boundary
my $h = pdl ([240, # blue -> green -> yellow ->
0]); # -> red
my $l = pdl ([0.6, 0.6]);
my $s = pdl ([0.8, 0.8]);
plscmap1n (256);
plscmap1l (0, $i, $h, $l, $s, pdl ([]));
}
my ($xm, $ym, $xM, $yM);
my $randn = 0;
my $rosen = 0;
sub main {
my $clev;
# Options data structure definition. */
my $pts = 500;
my $xp = 25;
my $yp = 20;
my $nl = 16;
my $knn_order = 20;
my $threshold = 1.001;
my $wmin = -1e3;
GetOptions ("npts=i" => \$pts,
"randn" => \$randn,
"rosen" => \$rosen,
"nx=i" => \$xp,
"ny=i" => \$yp,
"nlevel=i" => \$nl,
"knn_order=i" => \$knn_order,
"threshold=f" => \$threshold,
"help" => \$help);
if ($help) {
print (<<EOT);
$0 options:
--npts points Specify number of random points to generate [500]
--randn Normal instead of uniform sampling -- the effective
number of points will be smaller than the specified.
--rosen Generate points from the Rosenbrock function.
--nx points Specify grid x dimension [25]
--ny points Specify grid y dimension [20]
--nlevel Specify number of contour levels [15]
--knn_order order Specify the number of neighbors [20]
--threshold float Specify what a thin triangle is [1. < [1.001] < 2.]
EOT
push (@ARGV, "-h");
}
unshift (@ARGV, $0);
@title = ("Cubic Spline Approximation",
"Delaunay Linear Interpolation",
"Natural Neighbors Interpolation",
"KNN Inv. Distance Weighted",
"3NN Linear Interpolation",
"4NN Around Inv. Dist. Weighted");
@opt = (0., 0., 0., 0., 0., 0.);
$xm = $ym = -0.2;
$xM = $yM = 0.6;
plParseOpts (\@ARGV, PL_PARSE_PARTIAL);
$opt[2] = $wmin;
$opt[3] = $knn_order;
$opt[4] = $threshold;
# Initialize plplot
plinit ();
plseed (5489);
my ($x, $y, $z) = create_data ($pts); # the sampled data
my $zmin = min ($z);
my $zmax = max ($z);
my ($xg, $yg) = create_grid ($xp, $yp); # grid the data at
# the output grided data
plcol0 (1);
plenv ($xm, $xM, $ym, $yM, 2, 0);
plcol0 (15);
pllab ("X", "Y", "The original data sampling");
plcol0 (2);
plpoin ($x, $y, 5);
pladv (0);
plssub (3, 2);
for (my $k = 0; $k < 2; $k++) {
pladv (0);
for (my $alg = 1; $alg < 7; $alg++) {
$zg = plgriddata ($x, $y, $z, $xg, $yg, $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) {
for (my $i = 0; $i < $xp; $i++) {
for (my $j = 0; $j < $yp; $j++) {
if (not isfinite ($zg->slice ("$i,$j"))) {
# average (IDW) over the 8 neighbors
$zg->slice ("$i,$j") .= 0;
my $dist = 0;
for (my $ii = $i - 1; $ii <= $i + 1 && $ii < $xp; $ii++) {
for (my $jj = $j - 1; $jj <= $j + 1 && $jj < $yp; $jj++) {
my $zgij = $zg->slice ("$ii,$jj");
if ($ii >= 0 && $jj >= 0 && isfinite ($zgij)) {
$d = (abs ($ii - $i) + abs ($jj - $j)) == 1 ? 1. : 1.4142;
$zg->slice ("$i,$j") .= $zg->slice ("$i,$j")
+ $zgij / ($d * $d);
$dist += $d;
}
}
}
if ($dist != 0.) {
$zg->slice ("$i,$j") .= $zg->slice ("$i,$j") / $dist;
} else {
$zg->slice ("$i,$j") .= $zmin;
}
}
}
}
}
my $lzM = max ($zg);
my $lzm = min ($zg);
$lzm = min pdl ([$lzm, $zmin]);
$lzM = max pdl ([$lzM, $zmax]);
$lzm = $lzm - 0.01;
$lzM = $lzM + 0.01;
plcol0 (1);
pladv ($alg);
if ($k == 0) {
$clev = $lzm + ($lzM - $lzm) * sequence ($nl)/ ($nl - 1);
plenv0 ($xm, $xM, $ym, $yM, 2, 0);
plcol0 (15);
pllab ("X", "Y", $title [$alg - 1]);
plshades ($zg, $xm, $xM, $ym, $yM,
$clev, 1, 0, 1, 1, 0, 0, 0);
plcol0 (2);
} else {
$clev = $lzm + ($lzM - $lzm) * sequence ($nl)/ ($nl - 1);
cmap1_init ();
plvpor (0.0, 1.0, 0.0, 0.9);
plwind (-1.1, 0.75, -0.65, 1.20);
#
# For the comparition 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.0, -40.0);
plbox3 (0.0, 0, 0.0, 0, 0.5, 0,
"bntu", "X", "bntu", "Y", "bcdfntu", "Z");
plcol0 (15);
pllab ("", "", $title [$alg - 1]);
plot3dc ($xg, $yg, $zg, DRAW_LINEXY | MAG_COLOR | BASE_CONT, $clev);
}
}
}
plend ();
}
sub create_grid {
my ($px, $py) = @_;
my ($x, $y);
$x = $xm + ($xM - $xm) * sequence ($px) / ($px - 1);
$y = $ym + ($yM - $ym) * sequence ($py) / ($py - 1);
return ($x, $y);
}
sub create_data {
my $pts = shift;
my ($x, $y, $z);
# generate a vector of pseudo-random numbers 2*$pts long
my $t = zeroes(2*$pts);
plrandd($t); # use the PLplot standard random number generator for consistency with other demos
# Use every other random number for the x vect and y vect.
# This is done in this funny way to make the results identical
# to the c version, which calls plrandd once for x, then once for y in a loop.
my $xt = ($xM-$xm) * $t->mslice([0,(2*$pts)-2,2]);
my $yt = ($yM-$ym) * $t->mslice([1,(2*$pts)-1,2]);
if (not $randn) {
$x = $xt + $xm;
$y = $yt + $ym;
} else { # std=1, meaning that many points are outside the plot range
$x = sqrt(-2 * log($xt)) * cos(2 * pi * $yt) + $xm;
$y = sqrt(-2 * log($xt)) * sin(2 * pi * $yt) + $ym;
}
if (not $rosen) {
my $r = sqrt ($x ** 2 + $y ** 2);
$z = exp (- $r ** 2) * cos (2 * pi * $r);
} else {
$z = log ((1 - $x) ** 2) + 100 * ($y - $x ** 2) ** 2;
}
return ($x, $y, $z);
}
main ();