Friday, December 11, 2015

Perl for selecting data based on polygon


gnuplot > plot 'data.txt' u 1:2:3 pt 7 ps 2 t'', 'polygon.txt' u 1:2 w lines t'' 



gnuplot > plot 'inside_polygon.txt' u 1:2:3 pt 7 ps 2 t'', 'polygon.txt' u 1:2 w lines t''



data.txt
1 1 10
2 1 20
3 1 30
3 2 40
4 2 50
2 3 60
3 3 70
4 3 80
5 3 90
3 5 80
4 5 70
5 5 60
6 6 12
5 6 65
3 6 34
1 6 12

polygon.txt
1 3
3 1
6 3
4 6
2 5
1 3

#!/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 (FILE, 'data.txt');
while(  <FILE> ) {
chomp;
my @row = split;
push @angka, \@row;
}
close( FILE );
$nopoints=$#angka+1;
for my $i (0..$#angka)
{
push @px, $angka[$i][0];
push @py, $angka[$i][1];
push @pz, $angka[$i][2];
}

#### 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," > inside_polygon.txt" or die "$!\n";
for($i=0;$i <$#px+1;$i++)
{
$xxx[$i]=PtInPoly( \@poly, [$px[$i],$py[$i]]);
print OUT $px[$i],"\t";
print OUT $py[$i],"\t";
if ($xxx[$i] == 1)  
{
print OUT $pz[$i],"\n";
} else
{
print OUT 'NaN',"\n"; 
}
}

No comments: