Thanks, Matthew! This might be just what I need.
I hope your sins were not too grave...
I attach a script which uses your routine and plots a sinusoidal surface
both as generated and approximated by your routine. For this to work, you
need the newest version of PDL::Graphics::PLplot (v0.62) which I just put
on CPAN.
Regards,
Doug
[email protected]
Software Engineer
UCAR - COSMIC, Tel. (303) 497-2611
On Mon, 7 May 2012, Matthew Kenworthy wrote:
For my sins, I wrote a minimum curved surface interpolator based on
the same IDL routine. It tends to fall over for large sets of points,
but maybe it can help you! Let me know if this does the job.
Cheers,
Matt
On Mon, May 7, 2012 at 6:45 PM, Craig DeForest
<[email protected]> wrote:
I have done a few piecemeal things but there is a great need for irregular-grid
interpolation. It would be a great use for a Delaunay triangulator...
On May 7, 2012, at 10:43 AM, Doug Hunt wrote:
Hi all: I'm looking to draw a contour map in PDL/PLplot from scattered
soundings. So, before I use the PLplot function to draw a contour map from a
set of values on a fixed grid, I need to do some objective analysis, such as
Barnes interpolation.
Has anyone implemented such a thing in PDL? If not, are there any pointers to
software in other languages I could translate?
Thanks much!
--Doug Hunt
[email protected]
Software Engineer
UCAR - COSMIC, Tel. (303) 497-2611
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
--
Matthew Kenworthy / Assistant Professor / Leiden Observatory
Niels Bohrweg 2 (#463) / P.O. Box 9513 / 2300 RA Leiden / NL
#!/ops/tools/bin/perl
#
## Sample a 3D function, then grid the samples and do a contour plot
## of the regularized grid
#
use lib qw(/ops/tools/lib);
use strict;
use warnings;
use PDL;
use PDL::Graphics::PLplot;
my $n = 150;
my $nx = 31;
my $ny = 31;
my $ns = 15;
my $pi = 4*atan2(1,1);
my $x = random($n);
my $y = random($n);
my $z = 0.5 + 0.5 * (sin($x*2*$pi) * sin($y*2*$pi)) ** 3; # Bumps
my $xgrid = xvals($nx,$ny)/($nx-1);
my $ygrid = yvals($nx,$ny)/($ny-1);
my $zgrid = 0.5 + 0.5 * (sin($xgrid*2*$pi) * sin($ygrid*2*$pi)) ** 3; # Bumps
my $zinterp = min_curve_surf($x, $y, $z, $xgrid->flat, $ygrid->flat, 1)->reshape($xgrid->dims);
my $pl = PDL::Graphics::PLplot->new(DEV => 'pngcairo', FILE => 'contour2.png', PAGESIZE => [1000,800]);
$pl->shadeplot($zinterp, $ns, BOX => [$xgrid->minmax, $ygrid->minmax],
PALETTE => 'RAINBOW', PLOTTYPE => 'SHADE');
$pl->shadeplot($zgrid, $ns, COLOR => 'BLACK', PLOTTYPE => 'CONTOUR', CONTOURLABELS => 1);
$pl->close;
exit;
sub min_curve_surf {
#
# USAGE:
#
# ($zinterp) = min_curve_surf($x, $y, $z, $xgrid->flat, $ygrid->flat, $tps);
#
# (x, y, z) are 1D vectors defining the surface z = f(x,y)
# xgrid, ygrid are the points to sample to
# tcs sets whether to use TPS or MCS interpolation (see below)
# zinterp is a 1D returned vector
#
# EXAMPLE:
#
# $n = 15;
# $x = random(zeroes($n));
# $y = random(zeroes($n));
# $z = exp(-2*(($x-.5)^2 + ($y-.5)^2));
#
# $xgrid = xvals(10,10)/9; $ygrid = yvals(10,10)/9;
#
# ($zinterp) = min_curve_surf($x, $y, $z, $xgrid->flat, $ygrid->flat, 1);
# # THIN PLATE SPLINE = 1
# # MINIMUM CURVED SURFACE = 0
#
# print reshape $zinterp, dims($xgrid);
#
# REVISIONS:
#
# v04 - separated the back sub and decomp steps for clarification
# v03 - threaded lu_decomp2, gives 20% speed improvement
# v02 - is a working, confirmed version using unthreaded lu_decomp
# v01 - has debugging information so that it can be tested with IDL routines
#
use PDL::NiceSlice;
use PDL::MatrixOps;
my ($x,$y,$z, $xg, $yg, $tps) = @_;
my $n = nelem($x); # number of input points
my $tol = 1.0e-10;
# copying the IDL routine of the same name
# A minimum curvature spline surface is fitted to the data points
# described by X, Y, and Z. The basis function:
# C(x0,x1, y0,y1) = d^2 log(d^k),
# where d is the distance between (x0,y0), (x1,y1), is used,
# as described by Franke, R., "Smooth interpolation of scattered
# data by local thin plate splines," Computers Math With Applic.,
# v.8, no. 4, p. 273-281, 1982. k = 1 for minimum curvature surface,
# and 2 for Thin Plate Splines. For N data points, a system of N+3
# simultaneous equations are solved for the coefficients of the
# surface. For any interpolation point, the interpolated value
# is:
# F(x,y) = b(0) + b(1)*x + b(2)*y + Sum(a(i)*C(X(i),x,Y(i),y))
# The results obtained the thin plate spline (TPS) or the minimum
# curvature surface (MCS) methods are very similar. The only
# difference is in the basis functions: TPS uses d^2*alog(d^2),
# and MCS uses d^2*alog(d), where d is the distance from
# point (x(i),y(i)).
my $k = ($tps ? 1.0 : 0.5); # set the k exponent for tps
# make a square array that has n+3 sides
my $a = zeroes($n+3, $n+3);
# normalize z
my ($zmin, $zmax) = minmax($z);
my $zscale = 1.0 / ($zmax - $zmin);
# scale into 0-1 rect to enhance accuracy
my ($x0,$xmax) = minmax($x);
my ($y0,$ymax) = minmax($y);
my $xr = $xmax-$x0;my $yr =$ymax-$y0;
my $scale = 1. / ( ($xr>$yr) ? $xr : $yr);
my $xs = ($x - $x0) / $scale;
my $ys = ($y - $y0) / $scale;
### print "Fill in the distances between all points...";
# fill in the distances between all the points
my $tx1 = $xs * transpose(ones($n)) ;
my $tx2 = ones($n) * transpose($xs) ;
my $xmxn = $tx2 - $tx1;
my $ty1 = $ys * transpose(ones($n)) ;
my $ty2 = ones($n) * transpose($ys) ;
my $ymyn = $ty2 - $ty1;
my $d = ($xmxn*$xmxn) + ($ymyn*$ymyn);
### print "done\n";
$d += $tol*($tol > $d);
$d = $d * log($d) * $k;
$a(3:-1,0:$n-1) .= $d;
# add in other coefficients in the array
$a(0,0:$n-1) .= 1;
$a(1,0:$n-1) .= transpose($xs);
$a(2,0:$n-1) .= transpose($ys);
$a(3:-1,$n) .= 1;
$a(3:-1,$n+1) .= $xs;
$a(3:-1,$n+2) .= $ys;
my $b = zeroes($n+3);
$b(0:$n-1) .= ( ($z - $zmin) * $zscale); # sever z so that the ludecomp doesn't splash
### print "Start the lu decomposition...";
my ($lu, $perm, $parity) = lu_decomp2($a);
### print "done\n";
### print "Start the lu back substitution...";
my $c = lu_backsub($lu, $perm, $parity, $b);
### print "done\n";
# it's a whole square array, but just pull out the first line
$c = $c(,0;-);
# For min_curve_surf, divide c[3:*] by 2 cause we use d rather than
# the d^2 we use for TPS.
if (!$tps) {$c(3:-1) /= 2.0;}
# reapply scale factor to output grid
my $xpouts = ($xg - $x0) / $scale;
my $ypouts = ($yg - $y0) / $scale;
my $s = $c(0) + $c(1)*$xpouts + $c(2)*$ypouts;
### print "calcuating the final points...";
for (my $i =0; $i < $n; $i++) {
my $dd = pow($xpouts - $xs($i),2) + pow($ypouts - $ys($i),2);
$dd += $tol*($tol > $dd);
$s += ($dd * log($dd) * $c($i+3));
}
### print "done\n";
return(($s/$zscale) + $zmin);
}
1;
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
