Perl untuk memilih bagian dari peta 4D

Kode gnuplot di bawah ini adalah untuk mem-plot gambar 4D sebagaimana yang telah saya tunjukkan di sini. Tetapi untuk topik ini saya memanfaatkan select polygon untuk membatasi peta saya dengan kode Perl.

#!/bin/sh
gnuplot -persist << PLOT

set dgrid3d 300,300 gauss 1.5
dataFile='kontur.txt'

set table dataFile.'.grid'
splot dataFile u 1:2:3
unset table

set table dataFile.'.color'
splot dataFile u 1:2:4
unset table

set view 28,59
set hidden3d
set cbrange [1000:4000]
set palette defined (0 "red", 0.25 "orange", 0.5 "yellow", 0.75 "green", 1 "blue")
set cblabel "velocity (m/s)"
set autoscale cbfix
set pm3d

unset dgrid3d
set ticslevel 0
set xlabel 'X-COORD (m)' rotate parallel
set ylabel 'Y-COORD (m)' rotate parallel
set zlabel 'Depth (m)' rotate parallel
set zrange [400:800] reverse
set xrange[1:9]
set yrange[1:9]
splot sprintf('< paste %s.grid %s.color', dataFile, dataFile) u 1:2:3:7 with pm3d notitle

quit
PLOT

Kode Perl untuk memilih dalam polygon

#!/usr/bin/perl
use constant X         => 0;
use constant Y         => 1;
use constant PI        => atan2(0,-1);
use constant TWOPI    => 2*PI;

#### get data set
open (FILE2, 'kontur.txt.grid');
while( <FILE2> ) {
chomp;
my @s2 = split;

if ($_!~/#/)
{
if ($s2[i]>0)
{
push @px2, $s2[0];
push @py2, $s2[1];
push @pz2, $s2[2];
}
}
}
close( FILE2 );


#### printing
open OUT2," > kontur.grid" or die "$!\n";
for($i=0;$i<$#px2+1;$i++)
{
$neg2[$i]=$py2[$i+1]-$py2[$i];
if ($neg2[$i] < 0)  
{
print OUT2 $px2[$i],"\t";
print OUT2 $py2[$i],"\t";
print OUT2 $pz2[$i],"\t";
print OUT2 "\n";
} else
{
print OUT2 $px2[$i],"\t";
print OUT2 $py2[$i],"\t";
print OUT2 $pz2[$i],"\t";
}
print OUT2 "\n";
}


#### get data set
open (FILE, 'kontur.txt.color');
while( <FILE> ) {
chomp;
my @s = split;

if ($_!~/#/)
{
if ($s[i]>0)
{
push @px, $s[0];
push @py, $s[1];
push @pz, $s[2];
}
}
}
close( FILE );


#### get polygon
open (FILE1, 'polygon.txt');
while( <FILE1> ) {
chomp;
my @row1 = split;
push @angka1, \@row1;
}
close( FILE1 );
$nopoints1=$#angka1+1;
for my $i (0..$#angka1)
{
push @x, $angka1[$i][0];
push @y, $angka1[$i][1];
}


for($i=0;$i<$#x;$i++)
{
push(@{$poly[$i]}, $x[$i],$y[$i],); 
}


#### Detect Inside/Outside Polygon
sub mapAdjPairs (&@) {
    my $code = shift;
    map { local ($a, $b) = (shift, $_[0]); $code->() } 0 .. @_-2;
}

sub Angle{
    my ($x1, $y1, $x2, $y2) = @_;
    my $dtheta = atan2($y1, $x1) - atan2($y2, $x2);
    $dtheta -= TWOPI while $dtheta >   PI;
    $dtheta += TWOPI while $dtheta < - PI;
    return $dtheta;
}
sub PtInPoly{
    my ($poly, $pt) = @_;
    my $angle=0;

    mapAdjPairs{
        $angle += Angle(
            $a->[X] - $pt->[X],
            $a->[Y] - $pt->[Y],
            $b->[X] - $pt->[X],
            $b->[Y] - $pt->[Y]
        )
    } @$poly, $poly->[0];

    return !(abs($angle) < PI);
}
##################################


#### printing
open OUT," > kontur.color" or die "$!\n";
for($i=0;$i<$#px+1;$i++)
{
$xxx[$i]=PtInPoly( \@poly, [$px[$i],$py[$i]]);
$neg[$i]=$py[$i+1]-$py[$i];
if ($neg[$i] < 0)  
{
print OUT $px[$i],"\t";
print OUT $py[$i],"\t";
if ($xxx[$i] == 1)  
{
print OUT $pz[$i],"\n";
} else
{
print OUT 'NaN',"\n"; 
}
print OUT "\n";
} else

{
print OUT $px[$i],"\t";
print OUT $py[$i],"\t";
if ($xxx[$i] == 1)  
{
print OUT $pz[$i],"\n";
} else
{
print OUT 'NaN',"\n"; 
}
}


}


Kode gnuplot untuk menampilkan hasil:

#!/bin/sh
gnuplot -persist << PLOT
set view 28,59
set hidden3d
set cbrange [1000:4000]
set palette defined (0 "red", 0.25 "orange", 0.5 "yellow", 0.75 "green", 1 "blue")
set cblabel "velocity (m/s)"
set autoscale cbfix
set pm3d
set ticslevel 0
set xlabel 'X-COORD (m)' rotate parallel
set ylabel 'Y-COORD (m)' rotate parallel
set zlabel 'Depth (m)' rotate parallel
set zrange [400:800] reverse
set xrange[1:9]
set yrange[1:9]
splot sprintf('< paste kontur.grid kontur.color') u 1:2:3:6 with pm3d notitle

quit
PLOT