#!/usr/local/bin/perl
#
#dist.pl -- ３点、p0,p1,p2のそれぞれの間の距離をl0,l1,l2とするとき、
#			p1からl2への距離が小さいときにpqを省く。
#

require './glc-sub.pl';

$DEBUG = 0; # ログを出したいときは0以外に設定。
$T_LOG = 0; # トラックごとの圧縮結果を表示したいときは0以外に設定。

#更新履歴
#2007/09/22 作成開始。glc.plを流用してるので、いらないところもあるかも。
#2007/09/23 メインループの書き直し。すっきりしたような、してないような。
#2007/09/24 日付でスレッドを分割する機能を削除
#2007/09/29 移動距離が小さすぎるポイントも削除する機能を追加
#2007/10/21 時単位が変わったときは残すことに（ずっと直線でも１時間に１箇所ぐらいは残す）
#           ↑この機能は引き継いでない
#
#2008/04/16 １から作成しなおし
#2008/09/01 旅行に行って考えてみました。
#2008/09/20 圧縮ルーチン以外完成（ぉぃ）
#2008/09/23 いちおう完成？
#2008/09/28 つまらんバグ取り完了。
#2008/10/18 中身をすっきりさせて、１点ごとにソートするように変更。遅くなりすぎて非実用的。
#2008/10/19 とりあえず完了。
#2009/01/28 triangle.pl を流用して、面積から高さを求める方法で
#2009/01/30 点から線分への距離に変更

#カシミール３Ｄのpotファイル専用
#書式は以下
#DAT=WGS84/WGS84/0.000/0.00000000/0/0/0 
#ORD=POI/tex=/ido=034'41'24'5/kei=135'33'05'7/adr=1/alt=-5/col=12/pda=22-SEP-05/pti=21:44:54 
#ORD=POI/tex=/ido=034'41'26'3/kei=135'33'05'4/adr=1/alt=+9/col=12/pda=22-SEP-05/pti=21:45:23 
#ORD=POI/tex=/ido=034'41'25'8/kei=135'33'05'3/adr=2/alt=+9/col=12/pda=22-SEP-05/pti=21:47:20 
#ORD=POI/tex=/ido=034'41'25'8/kei=135'33'05'3/adr=2/alt=+11/col=12/pda=22-SEP-05/pti=21:47:24 
#                 緯度            経度        軌跡番号 高度  表示色    日付          時刻
#日付・時刻はUTC
#ORD=POI と tex= は固定にしておく

$t0=time();
$time0=localtime($t0);
print STDERR "$time0\n";

#許容誤差。これ以下なら省く。
#$MAX_ERROR = 4.36; #赤道での0.1秒角の斜辺長
#$MAX_ERROR = 4.330127019; #1辺10mの正三角形の高さ
$MAX_ERROR = 10;

$PI = 3.14159265358979;

$BIG_NUM = 40000000;

@tbl = ();
#@tbl=(POT,NO,LAT,LONG,TIME,ADTP);
#ハッシュとかほかに方法あるやろ
$POT  = 0;
$NO   = 1;
$PREV = 2;
$NEXT = 3;
$LAT  = 4;
$LONG = 5;
$TIME = 6;
$ADTP = 7; #advantage_of_delete_this_point
$DEL_F= 8;
$ADDR = 9;

#ヘッダー
$head = <>;
print $head;

#1行目
$i = 0;
$total_input  = 0;
$total_output = 0;

$first = <>;
chop $first;
$tbl[$i][$POT] = $first;
($ido ,$kei ,$address_1 ,$altitude ,$col ,$pda ,$pti) 
	= &get_variables($first);

$tbl[$i][$NO]   = 0;
$tbl[$i][$PREV] = 0;
$tbl[$i][$NEXT] = 1;
$tbl[$i][$LAT]  = (&dmspot2deg($ido)) / 180 * $PI; #lat
$tbl[$i][$LONG] = (&dmspot2deg($kei)) / 180 * $PI; #long
$tbl[$i][$TIME] = &pdapti2time($pda ,$pti); #time
$tbl[$i][$ADDR] = $address_1;

#2行目以降
while (<>){
#1点ログを読んで、
	chop ;
	$i++;
	$total_input++;

	$tbl[$i][$POT] = $_;
	($ido ,$kei ,$address ,$altitude ,$col ,$pda ,$pti) 
		= &get_variables($_);
	
	$tbl[$i][$NO]   = $i;
	$tbl[$i][$PREV] = $i-1;
	$tbl[$i][$NEXT] = $i+1;
	$tbl[$i][$LAT]  = &dmspot2deg($ido) / 180 * $PI; #lat
	$tbl[$i][$LONG] = &dmspot2deg($kei) / 180 * $PI; #long
	$tbl[$i][$TIME] = &pdapti2time($pda ,$pti); #time
	$tbl[$i][$ADDR] = $address;
	
#トラックの切れ目だったら
	if($address != $address_1){
		$trk_num++;
		print STDERR "Track No.$trk_num  " if($T_LOG);
		$tbl[$i-1][$NEXT] = $i-1;
		$stock = pop(@tbl);
		$num = @tbl;
		printf STDERR "track input %5d points , ",$num if($T_LOG);

#そのログを圧縮
		&compress_track();

#圧縮終わったら出力
		$num = &write_track();
		$total_output += $num;
		printf STDERR "output %5d points\n",$num if($T_LOG);
		
#前トラックで読みすぎたデータを次のトラックの頭に送る
		$address_1=$address;
		@tbl=();
		push (@tbl , $stock);
		$tbl[0][$NO]   = 0;
		$i=0;
	}
#ファイル全部読みこみ終わるまで繰り返し
}

#最後のトラックの処理
$i++;
$trk_num++;
$total_input++;

$num = @tbl;
print STDERR "Track No.$trk_num  " if($T_LOG);
printf STDERR "track input %5d points , ",$num if($T_LOG);

#最後のトラックを圧縮
&compress_track();

#圧縮終わったら出力
$num = &write_track();
$total_output += $num;
printf STDERR "output %5d points\n",$num if($T_LOG);

print STDERR "total input data $total_input point(s) , ";
print STDERR "total output data $total_output point(s)\n";

$t1=time();
$time1=localtime($t1);
$t=$t1-$t0;
print STDERR "$time1\n";
print STDERR "elapsed time $t\n";

exit;





#圧縮ルーチン
#&compress_track()
sub compress_track{
	
	if($#tbl < 2){return 0;}
	
	#各点を省いたときの誤差をセット
	&set_all_error();

	do{
		$del_num = 0;
		@idx = ();
		@idx = sort { @$a[$ADTP] <=> @$b[$ADTP] } @tbl;
			
		for ($j=0 ; $idx[$j][$ADTP] < $MAX_ERROR ; $j++){
			if($tbl[$idx[$j][$NO]][$ADTP] < $MAX_ERROR){
				#1点削除
				$del_num++;
				$tbl[$idx[$j][$NO]][$DEL_F] = 1;
				&del_point($idx[$j][$NO]);
			}
		}
		
	#最後に順番にならべなおす
	@tbl = sort { @$a[$NO] <=> @$b[$NO] } @tbl;#これいるの？
	
	}while ($del_num);
	
	return 1;
}
	
#出力
#$total_output += &write_log();
sub write_track{
	my($output_num) = 0;
	my($i) = 0;
	my($pot) = $tbl[$i][$POT];
	$pot =~ s/adr=\d*/adr=$trk_num/;
	print "$pot\n";
	for ( $i = 1 ; $i < @tbl ; $i++){
		if(!($tbl[$i][$DEL_F])){
			$output_num++;
			$pot = $tbl[$i][$POT];
			$pot =~ s/adr=\d*/adr=$trk_num/;
			print "$pot\n";
		}
	}
	return $output_num;
}

#debug用
sub write_table{
	$num    = @tbl;
	print STDERR "track input $num points \n" if ($DEBUG);
	$output_num= 0;
	for ( $i = 0 ; $i < @tbl ; $i++){
		print STDERR "$i,$tbl[$i][$NO],$tbl[$i][$PREV],$tbl[$i][$NEXT],$tbl[$i][$TIME],$tbl[$i][$ADTP],$tbl[$i][$DEL_F]\n" if ($DEBUG);
	}
	return $output_num;
}

#各点を省いたときの誤差をセット
#	&set_all_error();
sub	set_all_error{
	foreach $i (1..($#tbl-1)){
		&set_error($i);
	}

	$tbl[0][$ADTP] = $BIG_NUM;
	$tbl[$#tbl][$ADTP] = $BIG_NUM;

	return 1;
}

sub set_error{
	my($i)=@_;

	if($i == $#tbl){
		$tbl[$i][$ADTP] = $BIG_NUM;
		return 1;
	}
	elsif($i == 0){
		$tbl[$i][$ADTP] = $BIG_NUM;
		return 1;
	}

	my($p)=$tbl[$i][$PREV];
	my($n)=$tbl[$i][$NEXT];

	print STDERR "i=$i p=$p n=$n --- " if($DEBUG);
		$lat0 = $tbl[$p][$LAT] ;
		$long0= $tbl[$p][$LONG];
		$lat1 = $tbl[$i][$LAT] ;
		$long1= $tbl[$i][$LONG];
		$lat2 = $tbl[$n][$LAT] ;
		$long2= $tbl[$n][$LONG];
			
		$l0=(&hubeny ($lat0, $long0, $lat1, $long1))[0];
		$l1=(&hubeny ($lat1, $long1, $lat2, $long2))[0];
		$l2=(&hubeny ($lat0, $long0, $lat2, $long2))[0];
		if($l2 == 0){
			$tbl[$i][$ADTP] = ($l0+$l1)/2;
		}else{
			$dist = &heron($l0,$l1,$l2)/$l2;
			if(($l2**2+$l0**2)<$l1){$dist = $l0;}
			if(($l2**2+$l1**2)<$l0){$dist = $l1;}
		
			$tbl[$i][$ADTP] = $dist;
		}
	print STDERR "$dist\n" if($DEBUG);
	return 1;
}

#
sub del_point{
	my($i )=@_;
	my($p )=$tbl[$i][$PREV];
	my($n )=$tbl[$i][$NEXT];
	print STDERR "tipn $t,$i,$p,$n \n" if ($DEBUG);

	$tbl[$p][$NEXT]=$n;
	$tbl[$n][$PREV]=$p;
	
	&set_error($p);
	&set_error($n);
	$tbl[$i][$ADTP] = $BIG_NUM;
	
	print STDERR "del_point naino\n" if ($DEBUG);
	&write_table if ($DEBUG);
	
	return 1;
}

sub heron{
	my($a,$b,$c)=@_;
	my($s)=($a+$b+$c)/2;
	my($area)=($s*($s-$a)*($s-$b)*($s-$c));
	if($area>0){$area = sqrt($area);}else{$area = 0;}
	return $area;
}

# end of file