#!/usr/bin/perl
#
# geiger-map - map background radiation levels collected from gps-geiger
# Copyright (c) 2004 by Ian Kluft
# released under the GNU General Public License Version 2
#    http://www.kluft.com/~ikluft/opensource/GPLv2.txt
#
# some code derived from Rideplot, Copyright (c) 2000-2001 by Ian Kluft
# with contributions by others
#   http://rideplot.sourceforge.net/
#

use GD;                 # Graphics library
use POSIX qw( atan asin cos sin sqrt );
use Date::Calc qw( Delta_DHMS Decode_Date_US Decode_Month );
use Getopt::Long;

use strict;
use warnings;
use vars qw( $pi $deg $km_per_mile $debug $maxwidth $maxheight $brush_width
	@data $dist_m $maxlat $minlat $maxlon $minlon $max_x $max_y
	$maxgeiger $mingeiger $use_color @sp_color @sp_brush $image
	$img_x $img_y $white $black $brush_color $output_file
	%pts_by_lat %pts_by_lon $pixel_degrees_lat $pixel_degrees_lon
	@lat_keys %lat_keys @lon_keys %lon_keys $do_avg @color_gradient );

# convert degrees/min/sec to decimal degrees.
sub dms2deg
{
        my $str = shift;

        $debug and print STDERR "dms2deg coordinate: $str\n";

	if ( $str =~ /^-{0,1}[0-9]+\.[0-9]+$deg*$/o ) {
		return 0.0 + $str;
	}
	if ( $str =~ /^(-{0,1})([0-9]+)$deg*\s*([0-9]+\.[0-9]+)\'*$/o ) {
		return ( 0.0 + $2 + $3/60.0) * (($1 eq "-") ? -1.0 : 1.0 );
	}
	if ( $str =~ /^(-{0,1})([0-9]+)$deg*\s*([0-9]+)\'*\s*([0-9]+\.[0-9]+)\"*$/o ) {
		return ( 0.0 + $2 + $3/60.0 + $4/3600.0 ) * (($1 eq "-") ? -1.0 : 1.0 );
	}
	$debug and print STDERR "dms2deg: misformatted coordinate: $str\n";
	return undef;
}

# convert latitude/longitude to x/y
sub ll2xy
{
        my ( $lat, $lon, $maxlat, $minlat, $maxlon, $minlon, $max_x, $max_y )
                = @_;
        my ( $x, $y );
        $x = int (( $lon - $minlon ) / ( $maxlon - $minlon ) * $max_x )+$brush_width+1;
        $y = $max_y - ( int (( $lat - $minlat ) / ( $maxlat - $minlat ) * $max_y ))+$brush_width+1;
        $debug and print STDERR "ll2xy: lat=$lat lon=$lon max_x=$max_x max_y=$max_y x=$x y=$y\n";
        return ( $x, $y );
}

# convert degrees to radians for trig functions
sub deg2rad
{
        $_[0] * $pi / 180.0;
}

# convert radians to degrees for trig functions
sub rad2deg
{
        $_[0] * 180 / $pi;
}

# compute distance
sub dist
{
        my ( $lat1, $lon1, $lat2, $lon2 ) = @_;
        $debug and print STDERR "debug: dist ( $lat1, $lon1, $lat2, $lon2 )\n";

        # note: this great-circle forula should not be used near the poles
        # see http://users.netonecom.net/~rburtch/geodesy/datm_faq.html#2.1

        # Earth radius is estimated for each leg (Why? Because we can!)
        my $earth_radius = 6378.0 - 21.0 * sin(deg2rad(($lat1+$lat2)/2.0));

        # differences in lattitude and longitude
        my $dlat = abs($lat1-$lat2);
        my $dlon = abs($lon1-$lon2);

        # compute great-circle distance
        my $a = sin(deg2rad($dlat/2.0)) * sin(deg2rad($dlat/2.0))
                + cos(deg2rad($lat1))
                * cos(deg2rad($lat2))
                * sin(deg2rad($dlon/2.0)) * sin(deg2rad($dlon/2.0));
        my $c = 2.0 * asin(sqrt($a));

        # convert distance in radians on Earth-surface arc to meters
        $debug and print STDERR "distance (m) = "
                .($earth_radius * $c * 1000.0)
                .", c = $c, a = $a, dlat = $dlat, dlon = $dlon\n";
        return $earth_radius * $c * 1000.0;

}

# convert RGB hex string to red/green/blue integer tuples
sub rgb_hex2tuple
{
	my $hex = shift;

	if ( $hex =~ /^([0-9a-f]{2})([0-9a-f]{2})([0-9a-f]{2})$/i ) {
		return ( hex($1), hex($2), hex($3));
	} else {
		die "$0: $hex not in RGB hex string format\n";
	}
}

# initialize constants
$pi = 4.0 * atan(1);
$deg='°';
$km_per_mile = 1.609344;

# handle command-line parameters
$debug = 0;
$maxwidth = 800;
$maxheight = 600;
$brush_width=1;
$use_color = 1;
$do_avg = 0;

# initialize the graphing data array
@data = ();

# initialize distance accumulator
$dist_m = 0;

# command line
GetOptions ( "debug" => \$debug,
	"color:i" => \$use_color,               # Speed related colors. 2 inverts the colors.
	"output:s" => \$output_file,            # Output file name. STDOUT if not used
	"maxheight:i" => \$maxheight,           # Depreciated. For compatability reasons only.
	"maxwidth:i" => \$maxwidth,             # Depreciated. For compatability reasons only.
	"brush-width:s" => \$brush_width,
	"brush-color:s" => \$brush_color,
	"average|avg" => \$do_avg,
	"gradient|grad=s" => \@color_gradient, # series of rrggbb hex strings
);

if ( !@color_gradient ) {
	@color_gradient = ( "0000ff", "00ffff", "00ff00", "ffff00", "ff0000" );
}

my $datapoint = {};
while ( <> ) {
	chomp;
	if ( /^Time: ([0-9][0-9]:[0-9][0-9]:[0-9][0-9])/ ) {
		#print STDERR "Time: $1\n";
		$datapoint->{time} = $1;
		next;
	} elsif ( /^Position: ([0-9]+)\s+([0-9.]+)\s+([NS])\s+([0-9]+)\s+([0-9.]+)\s+([EW])/ ) {
		#print STDERR "Position: $1 $2 $3 $4 $5 $6\n";
		my $lat = dms2deg( $1." ".$2 );
		if ( $3 eq "S" ) {
			$lat = -$lat;
		}
		my $lon = dms2deg( $4." ".$5 );
		if ( $6 eq "W" ) {
			$lon = -$lon;
		}
		#print STDERR "Position: $lat $lon\n";
		$datapoint->{lat} = $lat;
		$datapoint->{lon} = $lon;

		# adjust high-water marks as needed
		if ( ! defined $maxlat ) {
			$maxlat = $minlat = $lat;
		} elsif ( $lat > $maxlat ) {
			$maxlat = $lat;
		} elsif ( $lat < $minlat ) {
			$minlat = $lat;
		}
		if ( ! defined $maxlon ) {
			$maxlon = $minlon = $lon;
		} elsif ( $lon > $maxlon ) {
			$maxlon = $lon;
		} elsif ( $lon < $minlon ) {
			$minlon = $lon;
		}

		next;
	} elsif ( /^(Geiger|Ur-hr): ([0-9.]+)/ ) {
		#print STDERR "Geiger: $1\n";
		my $geiger = $2;
		$datapoint->{geiger} = $geiger;

		# adjust high-water marks as needed
		if ( ! defined $maxgeiger ) {
			$maxgeiger = $mingeiger = $geiger;
		} elsif ( $geiger > $maxgeiger ) {
			$maxgeiger = $geiger;
		} elsif ( $geiger < $mingeiger ) {
			$mingeiger = $geiger;
		}

		next;
	} elsif ( ! /^\s*$/ ) {
		push ( @data, $datapoint );
		$datapoint = {};
	}
}

# create an index each for ordering data by latitude and longitude
# for searches for nearby data points
{
	my $i;
	for ( $i = 0; $i <= $#data; $i++ ) {
		my $lat = $data[$i]{lat} or next;
		my $lon = $data[$i]{lon} or next;
		if ( !defined $pts_by_lat{$lat}) {
			$pts_by_lat{$lat} = [];
		}
		if ( !defined $pts_by_lon{$lon}) {
			$pts_by_lon{$lon} = [];
		}
		push ( @{$pts_by_lat{$lat}}, $i );
		push ( @{$pts_by_lon{$lon}}, $i );
	}
	@lat_keys = sort keys %pts_by_lat;
	for ( $i = 0; $i <= $#lat_keys; $i++ ) {
		$lat_keys{$lat_keys[$i]} = $i;
	}
	@lon_keys = sort keys %pts_by_lon;
	for ( $i = 0; $i <= $#lon_keys; $i++ ) {
		$lon_keys{$lon_keys[$i]} = $i;
	}
}

# error checking
if ( $maxlat == $minlat ) {
        die "$0: no travel at all in latitude $maxlat - unable to scale\n";
}
if ( $maxlon == $minlon ) {
        die "$0: no travel at all in longitude $maxlon - unable to scale\n";
}

# generate the graph image object
# maintain the same aspect ratio as a real map of the track would

# compute size of image
{
	my $avg_lat = ( $maxlat + $minlat ) / 2;
	my $lat_scale = cos( $avg_lat * $pi / 180 );
	$max_x = int ($maxheight*(($maxlon-$minlon)/(($maxlat-$minlat)/$lat_scale)));
	$max_y = $maxheight;
	if ( $max_x > $maxwidth ) {
		$max_x = $maxwidth;
		$max_y = int($maxwidth*((($maxlat-$minlat)/$lat_scale))/($maxlon-$minlon));
	}
	$pixel_degrees_lat = ( $maxlat - $minlat ) / $max_y;
	$pixel_degrees_lon = ( $maxlon - $minlon ) / $max_x;
}

# allocate the image and colors
$img_x = ($max_x+$brush_width*2+2);
$img_y = ($max_y+$brush_width*2+2);
$image = new GD::Image($img_x, $img_y);
$image->setThickness($brush_width);

$white = $image->colorAllocate(255,255,255); # background
$black = $image->colorAllocate(0,0,0);
$image->transparent($white); # Make white the transparent color

if ( $use_color ) {
        my ( $i ); 
	my $colors = 100.0;
	my $grad_segs = scalar(@color_gradient);
        for ( $i=0; $i <= $colors; $i++ ) {
                my @rgb;
		#if ( $i > $colors/2 ) {
		#	@rgb = (255, int(255*(($colors-$i)/($colors/2))), 0);
		#} else {
		#	@rgb = (int(255*($i/($colors/2))), 255, 0);
		#}
                #push ( @sp_color, $image->colorAllocate(@rgb));

		if ( $i == 0 ) {
			@rgb = rgb_hex2tuple($color_gradient[0]);
		} elsif ( $i == $colors ) {
			@rgb = rgb_hex2tuple($color_gradient[$#color_gradient]);
		} else {
			my $seg = int( $i / $colors * ($grad_segs-1));
			my $seg_pos = ( $i - int($seg * ($colors/($grad_segs))))
				/ ($colors/($grad_segs));
			print STDERR "colors: i=$i seg=$seg pos=$seg_pos\n";
			my @seg_start = rgb_hex2tuple( $color_gradient[$seg]);
			my @seg_end = rgb_hex2tuple( $color_gradient[$seg+1]);
			my $red = int(($seg_end[0]-$seg_start[0]) * $seg_pos
				+ $seg_start[0]);
			my $green = int(($seg_end[1]-$seg_start[1]) * $seg_pos
				+ $seg_start[1]);
			my $blue = int(($seg_end[2]-$seg_start[2]) * $seg_pos
				+ $seg_start[2]);
			print STDERR "r=$red g=$green b=$blue\n";
			@rgb = ( $red, $green, $blue );
		}
                push ( @sp_color, $image->colorAllocate(@rgb));

                # make a brush for this color
                $brush_width == 1 and next;
                my $c_brush = new GD::Image($brush_width,$brush_width);
                my ( $c_fg, $c_bg );
                $c_fg = $c_brush->colorAllocate(@rgb);
                $c_bg = $c_brush->colorAllocate(255,255,255);
                $c_brush->transparent($c_bg);

                $c_brush->arc($brush_width,$brush_width,$brush_width/2,$brush_width/2,
                        0, 360, $c_fg );
                $c_brush->fill($brush_width/2,$brush_width/2,$c_fg);
                push ( @sp_brush, $c_brush );
        }
}

#
# iterate through all the data plotting points into the image
#
my $point;
my %pixel;
my ( $x, $y );
$debug and print STDERR "debug: brush width = $brush_width colors=$#sp_color brushes=$#sp_brush\n";

# first pass: loop through all the data points and save the data per pixel
for ( $point = 1; $point <= $#data; $point++ ) {
	( defined $data[$point]{geiger} ) or next;
	( defined $data[$point]{lat} ) or next;
	( defined $data[$point]{lon} ) or next;

	# save info for this point into the stack for its destination pixel
	( $x, $y ) = ll2xy($data[$point]{lat}, $data[$point]{lon},
		$maxlat, $minlat, $maxlon, $minlon, $max_x, $max_y );
	$debug and print STDERR "data: ($x,$y) <- ".$data[$point]{geiger}."\n";
	if ( !defined $pixel{$x}) {
		$pixel{$x} = {};
	}
	if ( !defined $pixel{$x}{$y}) {
		$pixel{$x}{$y} = [];
	}
	push ( @{$pixel{$x}{$y}}, $data[$point]{geiger});
}

# second pass: loop through all the pixels and generate the image
# note: keys are handled in shuffled order to mix overlapping pixels.
# (Pixels only overlap if brush size > 1.)
foreach $x ( sort { int(rand(3))-1 } keys %pixel ) {
	foreach $y ( sort { int(rand(3))-1 } keys %{$pixel{$x}} ) {
		if ( ! $use_color ) {
			$image->setPixel($x, $y, $black );
			$debug and print STDERR "plot ($x,$y) - black\n";
			next;
		}
		my ( $value, $seg_index );
		if ( $do_avg ) {
			# take average of data points for this pixel
			my ( $sum, $count, $i );
			$sum = 0.0;
			$count = 0;
			for ( $i = 0; $i < scalar @{$pixel{$x}{$y}}; $i++ ) {
				$sum += $pixel{$x}{$y}[$i];
				$count++;
			}
			$value = $sum/$count;
		} else {
			# take maximum of data points for this pixel
			my ( $i );
			$value = $pixel{$x}{$y}[0];
			for ( $i = 1; $i < scalar @{$pixel{$x}{$y}}; $i++ ) {
				if ( $pixel{$x}{$y}[$i] > $value ) {
					$value = $pixel{$x}{$y}[$i];
				}
			}
		}
		$seg_index = int( $value / $maxgeiger * $#sp_color );
		$debug and print STDERR "plot ($x,$y) - $value -> $seg_index\n";
		if ( $brush_width > 1 ) {
			$image->setBrush($sp_brush[$seg_index]);
			$image->setPixel($x, $y, gdBrushed );
		} else {
			my $seg_color = $sp_color[$seg_index];
			$image->setPixel($x, $y, $seg_color );
		}
	}
}

#  Open the track file if it's defined. 
if ( defined $output_file ) {
        if ( !open ( trk_img, ">$output_file" )) {
                print STDERR "$0: unable to open Track output file "
                        ."$output_file for writing: $!\n";
        }
	print trk_img $image->png;
	close trk_img ;

}