#!/usr/bin/perl

# Copyright © 2013 Guilhem Moulin <guilhem@fripost.org>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

use 5.010_000;
use warnings;
use strict;

our $VERSION = "0.1, August 9, 2013";

use Getopt::Long qw/:config posix_default no_ignore_case gnu_compat
                            bundling auto_version auto_help/;
use Pod::Usage;
use File::Temp;
use File::Basename;
use File::Spec;
use Inline 'C';

# http://geoinformatics.aalto.fi/doc/Geo-GDAL-1.9/html/
use Geo::GDAL;
use Geo::OSR;

# http://www.imagemagick.org/script/perl-magick.php
use Image::Magick;

my %config = ( progress   => Geo::GDAL->can('TermProgress_nocb') ?
                                sub { Geo::GDAL::TermProgress_nocb(@_[0,1]) } :
                                *Geo::GDAL::TermProgress
             , debug      => 0
             , rename     => 's/\.[^.]+$//'
             , zoom       => '1024x1024'
             , resampling => &resamplings('near')
             , border     => 10
             , bordercolor=> 'red'
             , fontsize   => 25
             , font       => 'Helvetica'
);


GetOptions( "q|quiet"  => sub { delete $config{progress} }
          , "d|debug"  => \$config{debug}
          , "rename=s" => sub { pod2usage(2) unless $_[1] =~ /^s(.).+\1[ixgcadlu]*$/;
                                $config{rename} = $_[1]
                              }
          , "s|size=s" => sub { if ($_[1] =~ /^(\d*\.?\d+)(%?)$/) {
                                    $config{zoom} = $2 ? 100. / $1 : $1;
                                }
                                elsif ($_[1] =~ /^\d+x\d+$/) {
                                    $config{zoom} = $_[1];
                                }
                                else {
                                    pod2usage(2);
                                }
                              }
          , "r|resampling=s"=> sub { $config{resampling} = &resamplings($_[1]) }
          , "b|border=s"=> sub { if ($_[1] =~ /^(?<b>\d*)(?<bc>:.*)?$/) {
                                     $config{border} = $+{b} if $+{b};
                                     $config{bordercolor} = $1 if $+{bc} =~ /^:(.+)/;
                                 }
                                 else {
                                      pod2usage(2);
                                 }
                               }
          , "f|font=s"  => sub { if ($_[1] =~ /^(?<fs>\d*)(?<fc>:.*)?(?<f>:[^:]*)?$/) {
                                     $config{fontsize} = $+{fs} if $+{fs};
                                     $config{fontcolor} = $1 if $+{fc} =~ /^:(.+)/;
                                     $config{font} = $1 if $+{f} =~ /^:(.+)/;
                                 }
                                 else {
                                      pod2usage(2);
                                 }
                               }
          , "s_srs=s" => sub { $config{s_srs} = Geo::OSR::SpatialReference::->new();
                               $config{s_srs}->SetFromUserInput( $_[1] ) }
          , "t_srs=s" => sub { $config{t_srs} = Geo::OSR::SpatialReference::->new();
                               $config{t_srs}->SetFromUserInput( $_[1] ) }
          , "man"     => sub { pod2usage(-exitstatus => 0, -verbose => 2) }
          )
    or pod2usage(2);
pod2usage(2) unless $#ARGV > 0;

my $index = pop;
my @mapset = map {{ filename => $_ }} @ARGV;

&debug( "Using mapset:", map {"  - ".$_->{filename}} @mapset );
&debug( "Building index: $index" );


# Find the source and target SRS
# We also store the corders points and the affine transformation, to
# avoid to compute it again later on.
my $s_srs;
foreach my $map (@mapset) {
    my $name = $map->{filename};

    &debug ("Opening '$name'");
    my $ds = Geo::GDAL::Open($name, 'ReadOnly') or exit 1;
    $map->{wkt} = $ds->GetProjection() || $config{s_srs}->ExportToWkt();

    unless ($config{s_srs}) {
        die "Error: Missing projection for '$name'.\n" unless $map->{wkt};

        $map->{srs} = Geo::OSR::SpatialReference::->new( $map->{wkt} );
        if ($s_srs) {
            unless ($s_srs->IsSame($map->{srs})) {
                my $msg = "Projection for '$name' differs from the previous input tiles'.";
                die "Error: " .$msg. "\n" unless $config{t_srs};
                warn "WARN: " .$msg. "\n";
            }
        }
        else {
            $s_srs = $map->{srs};
        }
    }

    &debug( "Using source SRS: ". ($config{s_srs} // $map->{srs}) );
    my $tx = Geo::OSR::CoordinateTransformation::->new(
                 $config{s_srs} // $map->{srs},
                 $config{t_srs} // $config{s_srs} // $s_srs
             );
    die "Error: Could compute transformation between source and target SRS.\n"
        unless $tx;

    my @geo = $ds->GetGeoTransform();
    $map->{transform} = \@geo;
    $map->{corners} = [ &pix2proj(\@geo,
                            ( [ 0,                    0]
                            , [ 0,                    $ds->{RasterYSize}-1 ]
                            , [ $ds->{RasterXSize}-1, 0 ]
                            , [ $ds->{RasterXSize}-1, $ds->{RasterYSize}-1 ]
                            ))
                      ];
    $tx->TransformPoints( $map->{corners} );
    die "Error: Could not transform point.\n"
        if grep !/^[-+]?\d*\.?\d+$/, (map @$_, @{$map->{corners}});


    { local $_ = basename $name; eval $config{rename}; $map->{text} = $_; }
}

sub pix2proj {
    my $geo = shift;
    map [ $geo->[0] + $_->[0]*$geo->[1] + $_->[1]*$geo->[2],
          $geo->[3] + $_->[0]*$geo->[4] + $_->[1]*$geo->[5] ],
        @_
}

$config{t_srs} //= $config{s_srs};
$config{s_wkt} = $config{s_srs}->ExportToWkt() if $config{s_srs};
$config{t_wkt} = $config{t_srs}->ExportToWkt();
&debug ("Using target SRS: $config{t_srs}");


# Compute the optimal zoom factor if the user is expecting a size.
if ($config{zoom} =~ /^(\d+)x(\d+)$/) {
    my ($xSize, $ySize) = ($1,$2);

    # Get the maximum resolution, and corners of the final tile;
    my ($xRes, $yRes) = (0,0);
    my ($ulx, $uly, $lrx, $lry);
    foreach (@mapset) {
        # TODO: What to do with multiple tiles?
        &debug ("Opening '$_->{filename}'");
        my $src_ds = Geo::GDAL::Open($_->{filename}, 'ReadOnly') or exit 1;
        my $dst_ds = Geo::GDAL::AutoCreateWarpedVRT( $src_ds
                                                   , $config{s_wkt} // $_->{wkt}
                                                   , $config{t_wkt} );
        $_->{transform2} = [ $dst_ds->GetGeoTransform() ];
        $config{ratio} = (0. + $dst_ds->{RasterXSize}) / $dst_ds->{RasterYSize};
        &debug( "Using ratio: $config{ratio}" );
        undef $src_ds;
        undef $dst_ds;

        $xRes = $_->{transform}->[1] if abs $xRes < abs $_->{transform}->[1];
        $yRes = $_->{transform}->[5] if abs $yRes < abs $_->{transform}->[5];

        $ulx = $_->{corners}->[0]->[0]
            unless defined $ulx and $ulx < $_->{corners}->[0]->[0];
        $uly = $_->{corners}->[0]->[1]
            unless defined $uly and $uly > $_->{corners}->[0]->[1];
        $lrx = $_->{corners}->[3]->[0]
            unless defined $lrx and $lrx > $_->{corners}->[3]->[0];
        $lry = $_->{corners}->[3]->[1]
            unless defined $lry and $lry < $_->{corners}->[3]->[1];
    }

    my $xRatio = (0.+$lrx-$ulx) / ($xRes*$xSize);
    my $yRatio = (0.+$lry-$uly) / ($yRes*$ySize);

    $config{zoom} = $xRatio < $yRatio ? $yRatio : $xRatio;
    map { $_->{transform2}->[1] = $_->{transform}->[1] * $config{zoom}
        ; $_->{transform2}->[5] = $_->{transform}->[5] * $config{zoom}
        } @mapset;
}
&debug( "Using zoom: $config{zoom}" );


# Working directory
my $tmpdir = File::Temp::tempdir( CLEANUP => 1 );
&debug ("Using working directory: '$tmpdir'.");

# Ensure the existence of the internal driver.
my $driver = 'GTiff';
my $driver_opts = { INTERLEAVE => 'BAND'
                  , TILED      => 'YES'
                  , COMPRESS   => 'NONE'
                  , BIGTIFF    => 'IF_NEEDED'
                  , PROFILE    => 'BASELINE'
                  , TFW        => 'YES'
                  };
my $driver_type = 'Byte';
die "Error: Unknown driver: '$driver'.\n"
    unless Geo::GDAL::GetDriverByName($driver);

foreach my $map (@mapset) {
    &debug ("Opening '$map->{filename}'");
    my $ds = Geo::GDAL::Open($map->{filename}, 'ReadOnly') or exit 1;
    &expand_rgba( \$ds, $map );
    my $f = &reproject( \$ds, $map );
    &annotate( $f, $map );
}


#######################################################################

# Print debug messages.
sub debug {
    print STDERR (join '', map {"DEBUG: " .$_. "\n"} @_) if $config{debug};
}

# Create a temporary file. The file is created to defeat race
# conditions, but the file handle is closed right away.
sub tempfile {
    my ($fh, $name) = File::Temp::tempfile(DIR => $tmpdir, UNLINK => 1, @_);
    close $fh;
    &debug( "Creating file: '$name'" );
    return $name;
}

# The available resampling methods. The strings can be found
# with keys %Geo::GDAL::RESAMPLING_STRING2INT;
sub resamplings {
    my $r = shift;
    my %h = ( near => 'NearestNeighbour'
            , bilinear => 'Bilinear'
            , cubic => 'Cubic'
            , cubicspline => 'CubicSpline'
    );
    $h{$r};
}

# Expand a dataset to RGBA. This code is mostly a perl translation of pct2rgb(1).
sub expand_rgba {
    my $ds = shift;
    my $map = shift;
    my $src_ds = $$ds;

    return unless $src_ds->{RasterCount} == 1;

    my $src_band = $src_ds->GetRasterBand(1);
    my $ct = $src_band->GetRasterColorTable() or return;

    &debug("Expansion to RGBA required.");

    my $dst_nb = 4;     # RGBA
    my $src_nc = 256;   # That's the best we can do, as we're using bytes

    my $lookup = [ map '', (0 .. $src_nc-1) ];
    $src_nc = $ct->GetCount() if $ct->GetCount() < $src_nc;
    for (my $i = 0; $i < $src_nc; $i++) {
        my @e = $ct->GetColorEntryAsRGB( $i );
        $lookup->[$_] .= chr $e[$_] for (0 .. $dst_nb-1);
    }

    my $dst_filename = &tempfile(SUFFIX => '.tif');
    my ($xSize, $ySize) = $src_ds->Size();
    my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $dst_filename
                                                            , $xSize
                                                            , $ySize
                                                            , $dst_nb
                                                            , $driver_type
                                                            , $driver_opts
                                                            );

    $dst_ds->SetProjection( $map->{wkt} );
    $dst_ds->SetGeoTransform( $map->{transform} );
    $dst_ds->SetGCPs( $src_ds->GetGCPs(), $src_ds->GetGCPProjection() )
        if $src_ds->GetGCPCount();

    &{$config{progress}}(0.0) if $config{progress};
    for (my $y = 0; $y < $ySize; $y++) {
        my $src_data = $src_band->ReadRaster(0, $y, $xSize, 1);

        for (my $b = 0; $b < $dst_nb; $b++) {
            local $_ = $src_data;
            # We use an inlined C function to speed up things a bit.
            # Essentially what we're doing is a splice through the
            # lookup table.
            # Using strings is much more efficients that huge perl
            # array of scalars. A pure perl version, not so much slower
            # than our C code, can be obtained using transliteration:
            #   my $src_str = join '', map chr, (0 .. $src_nc-1);
            #   eval "tr/\Q$src_str\E/\Q$lookup->[$b]\E/";
            &subst ( $_, $lookup->[$b], $xSize );
            $dst_ds->Band(1+$b)->WriteRaster(0, $y, $xSize, 1, $_);
        }

        &{$config{progress}}( (1+$y) / $ySize )
            if $config{progress};
    }

    # Erase (hence flush and close) the previous dataset
    undef $src_ds;
    $$ds = $dst_ds;
    return $dst_filename;
}


# Reproject the map. This code is inspired from
# http://jgomezdans.github.io/gdal_notes/reprojection.html
sub reproject {
    my $ds = shift;
    my $map = shift;
    my $src_ds = $$ds;

    my @geo = @{$map->{transform}};
    my ($ul,$lr) = @{$map->{corners}}[0,3];

    # TODO: it might be xSize that's determined from ySize
    my $dst_xSize = int( ($lr->[0]-$ul->[0]) / $map->{transform2}->[1] );
    my $dst_ySize = int( $dst_xSize / $config{ratio} );

    my ($src_xSize, $src_ySize) = $src_ds->Size();
    &debug("Reducing map from ${src_xSize}x${src_ySize} to ${dst_xSize}x${dst_ySize}.");

    my $dst_filename = &tempfile(SUFFIX => '.tif');
    my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $dst_filename
                                                            , $dst_xSize
                                                            , $dst_ySize
                                                            , $src_ds->{RasterCount}
                                                            , $driver_type
                                                            , $driver_opts
                                                            );
    $dst_ds->SetProjection( $config{t_wkt} );
    $dst_ds->SetGeoTransform( $map->{transform2} );

    my $res = Geo::GDAL::ReprojectImage( $src_ds, $dst_ds,
                                       , $map->{wkt}, $config{t_wkt}
                                       , $config{resampling}, undef, undef
                                       , $config{progress} );
    die "Error: Failed to reproject map " . join (', ', map {"'".$_."'"}
                                                        $src_ds->GetFileList() )
                                          . ".\n"
        if $res;

    # Erase (hence flush and close) the previous dataset
    undef $src_ds;
    return $dst_filename;
}


# Annotate the map (add border & caption).
sub annotate {
    my $src_filename = shift;
    my $map = shift;

    my $img = Image::Magick::->new();

    $img->Read($src_filename);
    my ($xSize, $ySize) = $img->Get(qw/width height/);

    &debug( "Border: $config{border}:$config{bordercolor}" );
    &debug( "Font: $config{fontsize}:". ($config{fontcolor} // $config{bordercolor}).
            ":$config{font}" );
    &debug( "Text: $map->{text}");
    # To get the list of fonts:
    #   &debug( "Supported fonts:", map {"\t".$_} $img->QueryFont() );

    my @draw = ( fill => 'none'
               , stroke => $config{bordercolor}
               , strokewidth => $config{border}
               );
    my @ann = ( fill => $config{fontcolor} // $config{bordercolor}
              , strokewidth => 0
              , font => $config{font}
              , pointsize => $config{fontsize}
              , gravity => 'Center'
              , text => $map->{text}
              );
    my @points;

    if ($config{t_srs}->IsSame($config{s_srs} // $map->{srs})) {
        push @draw, primitive => 'rectangle';
        push @points, [ 0, 0 ], [ $xSize-1, $ySize-1 ];
    }
    else {
        # The projection changed. Since the map might have been rotated,
        # we draw a polygon instead of a rectangle. Also, instead of
        # adding the border *then* reproject we use a mask in order to
        # have nice anti-aliasing and better looking fonts.

        my ($ul1,$ll1,$ur1,$lr1) = @{$map->{corners}};
        my ($rx,$ry) = @{$map->{transform2}}[1,5];
        my ($ul2,$ll2,$ur2,$lr2) = &pix2proj( $map->{transform2},
                                              ( [ 0,        0]
                                              , [ 0,        $ySize-1 ]
                                              , [ $xSize-1, 0 ]
                                              , [ $xSize-1, $ySize-1 ]
                                            ));
        push @draw, primitive => 'polygon';
        push @points, [ int(($ul1->[0]-$ul2->[0]) / $rx)
                      , int(($ul1->[1]-$ul2->[1]) / $ry) ]
                    , [ int(($ll1->[0]-$ll2->[0]) / $rx)
                      , int(($ll1->[1]-$ll2->[1]) / $ry) + $ySize ]
                    , [ int(($lr1->[0]-$lr2->[0]) / $rx) + $xSize
                      , int(($lr1->[1]-$lr2->[1]) / $ry) + $ySize ]
                    , [ int(($ur1->[0]-$ur2->[0]) / $rx) + $xSize
                      , int(($ur1->[1]-$ur2->[1]) / $ry) ]
                    ;
    }
    $img->Draw( @draw, points => join (' ', map {join ',', @$_} @points) );
    $img->Annotate( @ann );

    my $dst_filename = &tempfile(SUFFIX => '.tif');
    $img->Write($dst_filename);

    # The georeference didn't change when annotating. Hence we simply
    # link the former world file.
    my ($src_tfw, $dst_tfw) = map { my ($name,$path) = fileparse($_, qr/\.tiff?$/);
                                    File::Spec->catfile($path, $name.'.tfw');
                                  }
                                  ($src_filename, $dst_filename);

    if (-f $src_tfw and not -f $dst_tfw) {
        link $src_tfw, $dst_tfw;
        &debug( "Linking worldfiles: $src_tfw -> $dst_tfw" );
    }

    undef $img;
    $map->{thumbnail} = $dst_filename;
}


__END__

__C__

void subst (unsigned char *data, unsigned char *lookup, int length) {
    int i;
    for (i = 0; i < length; i++)
        data[i] = lookup[ data[i] ];
}