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:
Post a Comment