Friday, March 6, 2015

Boundary Detection Perl


#!/usr/bin/perl
$file=$ARGV[0];
my $pi = 3.14159265358979;
my @data;
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
push @data, \@row;
}
close( FILE );

for my $i (0..$#data)
{
push @pri, $data[$i][0];
push @sec, $data[$i][1];
}

for my $i (0..$#data-1)
{
$dx[$i]=$data[$i][1]-$data[$i+1][1];   #col 2 diff
if($dx[$i] != 0) 
{
push @ind,$i,$i+1;
}
}

push @ind1, 0,@ind, $#data;
for my $i (0..$#data)
{
if($ind1[$i] != 0) 
{
push @pria, $pri[$ind1[$i]];
push @seca, $sec[$ind1[$i]]
}
}

push @x1,$pri[0],@pria;
push @y1,$sec[0],@seca;

### try this option first
#$x0=0.1+(min(@pri)+max(@pri))/2;
#$y0=0.2+(min(@sec)+max(@sec))/2;

### if you see 'saw' shape, privide center point manually
$x0=800.1;
$y0=3530.2;

#calculate azimuth of each point relative to the center point
for my $i (0..$#x1)
{
$dy[$i]=abs($y1[$i]-$y0);
$dx[$i]=abs($x1[$i]-$x0);
$theta[$i]=rad_to_deg(atan($dy[$i]/$dx[$i]));
if ($x1[$i] < $x0 && $y1[$i]>$y0 )
{
$thetacor[$i]=180-$theta[$i];
}elsif ($x1[$i] < $x0 && $y1[$i]<$y0 )
{
$thetacor[$i]=180+$theta[$i];
}elsif ($x1[$i] > $x0 && $y1[$i]<$y0 )
{
$thetacor[$i]=360-$theta[$i];
}else
{
$thetacor[$i]=$theta[$i];
}
push(@{$dataku[$i]}, $x1[$i],$y1[$i],$thetacor[$i]);    #push multidimensional array
}

open OUT," > out.txt" or die "$!\n";
my @sorted = sort { $a->[2] <=> $b->[2] } @dataku;  ##sort by col 2 (azimuth)

#Print output
for my $i (0..$#sorted)
{
for my $j (0..$#{$sorted[0]}-1)
{
   print OUT  "$sorted[$i][$j]","\t";
}
   print OUT  "\n";
}

for my $j (0..$#{$sorted[0]}-1)   #close polygon
{
print OUT  "$sorted[0][$j]","\t";
}

print OUT  "\n";


sub max {
    splice(@_, ($_[0] > $_[1]) ? 1 : 0, 1);
    return ($#_ == 0) ? $_[0] : max(@_);
}

sub min {
    splice(@_, ($_[0] > $_[1]) ? 0 : 1, 1);
    return ($#_ == 0) ? $_[0] : min(@_);
}

sub deg_to_rad { ($_[0]/180) * $pi }
sub rad_to_deg { ($_[0]/$pi) * 180 }
sub asin { atan2($_[0], sqrt(1 - $_[0] * $_[0])) }
sub acos { atan2( sqrt(1 - $_[0] * $_[0]), $_[0] ) }
sub tan  { sin($_[0]) / cos($_[0])  }
sub atan { atan2($_[0],1) };

########## PLOT#############################################
open my $GP, '|-', 'gnuplot -persist';
print {$GP}  <<'__GNUPLOT__';
set multiplot;                          # get into multiplot mode
set size 1,0.5; 
set origin 0.0,0.5;
set xrange[700:1150]
set yrange[3480:3570]
plot 'PRIM_SEC2.TXT' using 1:2 with points pointtype 0 pointsize 1 title ""

set origin 0.0,0.0; 
set xrange[700:1150]
set yrange[3480:3570]
plot 'out.txt' using 1:2 with lines title ""
__GNUPLOT__
close $GP;


Kode di bawah ini adalah versi 2:
#!/usr/bin/perl
$file=$ARGV[0];
my @data;
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
push @data, \@row;
}
close( FILE );

for my $i (0..$#data)
{
push @pri, $data[$i][0];
push @sec, $data[$i][1];
}

for my $i (0..$#data-1)
{
$dx[$i]=$data[$i][1]-$data[$i+1][1];   #col 2 diff
if($dx[$i] != 0) 
{
push @ind,$i,$i+1;
}
}

push @ind1, 0,@ind, $#data;
for my $i (0..$#data)
{
if($ind1[$i] != 0) 
{
push @pria, $pri[$ind1[$i]];
push @seca, $sec[$ind1[$i]]
}
}

push @x1,$pri[0],@pria;
push @y1,$sec[0],@seca;


for my $i (0..$#x1/2)
{
push @{$left[$i]},$x1[$i*2],$y1[$i*2];    #push multidimensional array
push @{$right[$i]},$x1[($i*2)+1],$y1[($i*2)+1];   #push multidimensional array
}

my @sortedleft= sort { $a->[1] <=> $b->[1] } @left;  ##sort accending by col 2
my @sortedright = sort { $b->[1] <=> $a->[1] } @right;  ##sort decending by col 2

push @sorted, @sortedleft, @sortedright;
open OUT," > out.txt" or die "$!\n";
#Print output
for my $i (0..$#sorted)
{
for my $j (0..$#{$sorted[0]})
{
   print OUT   "$sorted[$i][$j]","\t";
}
   print OUT  "\n";
}

for my $j (0..$#{$sorted[0]})   #close polygon
{
print OUT  "$sorted[0][$j]","\t";
}

print OUT  "\n";

########## PLOT#############################################
open my $GP, '|-', 'gnuplot -persist';
print {$GP}  <<'__GNUPLOT__';
set multiplot;                          # get into multiplot mode
set size 1,0.5; 
set origin 0.0,0.5;
set xrange[700:1150]
set yrange[3480:3570]
plot 'PRIM_SEC.TXT' using 1:2 with points pointtype 0 pointsize 1 title ""

set origin 0.0,0.0; 
set xrange[700:1150]
set yrange[3480:3570]
plot 'out.txt' using 1:2 with lines title ""
__GNUPLOT__
close $GP;



No comments: