#!/usr/bin/perl
#
# 数値地図(標高)からUSGS海岸線データを
# Nowral
# 98.12.12
#
@filename = @ARGV;
foreach $oldargv (@filename) {
open(ARGVIN, "<$oldargv");

$newargv = $oldargv . '.dat';
open(ARGVOUT, ">$newargv");
select(ARGVOUT);
&MacPerl'SetFileInfo('JEDT', 'TEXT', $newargv);

$_ = <ARGVIN>; #ヘッダー行
/^[\s\d]{23}(\d{3})(\d{3})(\d{7})(\d{7})(\d{7})(\d{7})/;
$nr   = $1;
$lat2 = &dms2deg($3);
$lon0 = &dms2deg($4);
$lat0 = &dms2deg($5);
$lon2 = &dms2deg($6);
$dlat = ($lat2-$lat0) / $nr;
$dlon = ($lon2-$lon0) / $nr;

#ビット列作り
#$slb[$y] < j番目の走査線 全データ格納
#2 壁 1 陸 0 海
$slb[0] = "2" x ($nr+2);
$y = 1;
while (<ARGVIN>) {
    if(/\d+/) { #空行飛ばし
        s/\d{4}(\d{5})//;   #先頭の図番号を削除
        for(; $y < $1; ++$y) { $slb[$y] = "2","0" x $nr,"2"; } #海行省略
        s/[\n\r]//g;
        s/([\-\d]\d{4})/0+($1>-9999)/eg; #境界 <>等高線
        $slb[$y] = "2" . $_ . "2";
        ++$y;
    }
}
$slb[$y] = "2" x ($nr+2);
#foreach (@slb) {print $_,"\n";}

#方向テーブル
#北から時計回り
#    0
#    |
# 3- @ -1
#    |
#    2
@dx = ( 0, 1, 0,-1);
@dy = (-1, 0, 1, 0);

#plaquetteテーブル
# +-+-+
# |3|0|
# +-@-+
# |2|1|
# +-+-+
@px = ( 1, 1, 0, 0);
@py = ( 0, 1, 1, 0); 

#発端
$sc = 0;

#上辺
$x = 0;
$y = 0;
$pb = substr($slb[$y+$py[1]],$x+$px[1],1);
for($x = 1 ; $x<$nr; ++$x) {
    $nb = substr($slb[$y+$py[1]],$x+$px[1],1);
    if($pb==1 && $nb==0) {
        $spx[$sc] = $x;
        $spy[$sc] = $y;
        $spd[$sc] = 2;
        ++$sc;
    }
    $pb = $nb;
}

#右辺
$x = $nr;
$y = 0;
$pb = substr($slb[$y+$py[2]],$x+$px[2],1);
for($y = 1 ; $y<$nr; ++$y) {
    $nb = substr($slb[$y+$py[2]],$x+$px[2],1);
    if($pb==1 && $nb==0) {
        $spx[$sc] = $x;
        $spy[$sc] = $y;
        $spd[$sc] = 3;
        ++$sc;
    }
    $pb = $nb;
}

#下辺
$x = $nr;
$y = $nr;
$pb = substr($slb[$y+$py[3]],$x+$px[3],1);
for($x=$nr-1 ; $x>0; --$x) {
    $nb = substr($slb[$y+$py[3]],$x+$px[3],1);
    if($pb==1 && $nb==0) {
        $spx[$sc] = $x;
        $spy[$sc] = $y;
        $spd[$sc] = 0;
        ++$sc;
    }
    $pb = $nb;
}

#左辺
$x = 0;
$y = $nr;
$pb = substr($slb[$y+$py[0]],$x+$px[0],1);
for($y=$nr-1 ; $y>0; --$y) {
    $nb = substr($slb[$y+$py[0]],$x+$px[0],1);
    if($pb==1 && $nb==0) {
        $spx[$sc] = $x;
        $spy[$sc] = $y;
        $spd[$sc] = 1;
        ++$sc;
    }
    $pb = $nb;
}

#内側、閉曲線
for($y=1 ; $y<$nr; ++$y) {
    $x = 0;
    $pb = substr($slb[$y+$py[1]],$x+$px[1],1);
    for($x=1 ; $x<$nr; ++$x) {
        $nb = substr($slb[$y+$py[1]],$x+$px[1],1);
        if($pb==1 && $nb==0) {
            $spx[$sc] = $x;
            $spy[$sc] = $y;
            $spd[$sc] = 2;
            ++$sc;
        }
        $pb = $nb;
    }
}

#探索
$lcl = -1; #海岸線長 length of coast line

spl:foreach $sc (0..$#spx) { #出発点ループ
    
    $x = $spx[$sc];
    $y = $spy[$sc];
    if($x == 0 && $y == 0) { next spl; }
    
    ++$lcl;
    $clx[$lcl] = 0;
    $cly[$lcl] = 0;
    ++$lcl;
    $clx[$lcl] = $x;
    $cly[$lcl] = $y;
    $pd = $spd[$sc];
    $lcl0 = $lcl;

qcl:while(1) { #探索ループ
        $nd = ($pd+3) % 4;
        $rb = substr($slb[$y+$py[$nd]],$x+$px[$nd],1);
        if($rb == 2) { last qcl; } #壁にぶつかって終わり
        if($rb == 0) {
            $nd++;
            $nd %= 4;
            $rb = substr($slb[$y+$py[$nd]],$x+$px[$nd],1);
            if($rb == 0) {
                $nd++;
                $nd %= 4;
            }
        }
        if($nd != $pd) {
            ++$lcl;
            $clx[$lcl] = $x;
            $cly[$lcl] = $y;
            $pd = $nd;
            if($x==$spx[$sc] && $y==$spy[$sc]) { #曲線が閉じて終わり
                $x = $clx[$lcl0+1];
                $y = $cly[$lcl0+1];
                last qcl;
            }
        }
        
        $x += $dx[$pd];
        $y += $dy[$pd];

        foreach $i ($sc+1..$#spx) {
            if($x==$spx[$i] && $y==$spy[$i]) {
                $spx[$i] = 0;
                $spy[$i] = 0;
            }
        }
    }
    ++$lcl;
    $clx[$lcl] = $x;
    $cly[$lcl] = $y;

#近道
#中点処理
    $wc = 0;
    $x2 = $clx[$lcl0];
    $y2 = $cly[$lcl0];
    $x3 = $clx[$lcl0+1];
    $y3 = $cly[$lcl0+1];
    $d23 = sqrt(($x2-$x3)**2 + ($y2-$y3)**2);
    $tx3 = $x2;
    $ty3 = $y2;
    
    foreach $i ($lcl0+2 .. $lcl) {
        $x1 = $x2;
        $y1 = $y2;
        $x2 = $x3;
        $y2 = $y3;
        $x3 = $clx[$i];
        $y3 = $cly[$i];
        
        $d12 = $d23;
        $d23 = sqrt(($x2-$x3)**2 + ($y2-$y3)**2);
#        print $d12,",",$d23,"\n";

        $tx1 = $tx3;
        $ty1 = $ty3;
        if($d12==1 || $d23==1) {
            $tx2 = ($x2 + $x1) / 2;
            $ty2 = ($y2 + $y1) / 2;
            $tx3 = ($x2 + $x3) / 2;
            $ty3 = ($y2 + $y3) / 2;
        }
        else {
            $tx2 = $x2 + ($x1 - $x2) / $d12 / 2;
            $ty2 = $y2 + ($y1 - $y2) / $d12 / 2;
            $tx3 = $x2 + ($x3 - $x2) / $d23 / 2;
            $ty3 = $y2 + ($y3 - $y2) / $d23 / 2;
        }
        
        if($tx2 != $tx1 || $ty2 != $ty1) {
            $wx[$wc] = $tx1;
            $wy[$wc] = $ty1;
            ++$wc;
        }        
        $wx[$wc] = $tx2;
        $wy[$wc] = $ty2;
        ++$wc;
    }
     
    $wx[$wc] = $tx3;
    $wy[$wc] = $ty3;
    ++$wc;
    $wx[$wc] = $clx[$lcl];
    $wy[$wc] = $cly[$lcl];
    
    foreach $i (0..$wc) {
        $clx[$lcl0 + $i] = $wx[$i];
        $cly[$lcl0 + $i] = $wy[$i];
    }
    $lcl = $lcl0 + $wc;

#省略
    $wc = $lcl0 - 1;
    $px = $clx[$lcl0];
    $py = $cly[$lcl0];
    $pa = 0;
    foreach $i ($lcl0+1..$lcl) {
        $nx = $clx[$i];
        $ny = $cly[$i];
        if($nx != $px) {$na = ($ny - $py) / ($nx - $px); }
        else { $na = $nr + 1; }
        if($pa != $na) {
            ++$wc;
            $clx[$wc] = $px;
            $cly[$wc] = $py;
            $pa = $na;
        }
        $px = $nx;
        $py = $ny;
    }
    ++$wc;
    $clx[$wc] = $px;
    $cly[$wc] = $py;
    $lcl = $wc;
    
}

#出力　decode
foreach $i (0..$lcl) {
    $x = $clx[$i];
    $y = $cly[$i];
    if($x==0 && $y==0) {
        print "# -b\n";
    }
    else {
        $lon = $lon0 + $x*$dlon;
        $lat = $lat0 + $y*$dlat;
        printf "%.5f\t%.5f\n",$lon,$lat;
#        print $x,",",$y,"\n";
    }
}

}
&MacPerl'Quit(2);
die "The Unhappy End";


sub dms2deg {
    local($_) = @_;
    s|(\d\d\d)(\d\d)(\d\d)|$1 + $2/60 + $3/3600|e;
    $_;
}