#!/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/>.

our $VERSION = "0.3, August 13, 2013";

use 5.010_000;
use warnings;
use strict;

=head1 NAME

mkindex.pl - Merge a mapset into an index file

=head1 SYNOPSIS

B<mkindex.pl> [B<-q>] [B<-d>] [B<--rename> I<perlexpr>] [B<-s> I<geometry>]
[B<-r> I<method>] [B<-b> I<border>] [B<-f> I<font>] [B<--s_srs> I<SRS>]
[B<--t_srs> I<SRS>] I<map> [... I<map>] I<index>

B<mkindex.pl> B<--man>

=head1 DESCRIPTION

Merge all given I<map>s into an index file, and annotate this index with
tile borders and names.

If required, each map is first expanded into RGBA, to get a decent
resampling. This expansion can take a significant amount of space. An
alternative, much slower but which require nearly no extra space, is to
create expanded VRTs:

    for x in *.tif; do gdal_translate -expand rgba -of VRT $x ${x%.tif}.vrt; done
    mkindex.pl *.vrt index.png

=head1 OPTIONS

=over 8

=item B<-q>, B<--quiet>

Be quiet.

=item B<-d>, B<--debug>

Debug mode.

=item B<--rename>=I<perlexpr>

A Perl expression which is expected to modify the $_ string in Perl (set
to each filename), to what should be displayed on the final index. The
default is 's/\.[^.]+$//' that is, only the file extension is stripped
away. See the B<EXAMPLES> section for other examples.

=item B<-s> I<width>xI<height>, B<--size=>I<width>xI<height>

The size of the output index (default: 1024x0). If I<width> or I<height>
is 0, the other dimension will be guessed from the computed resolution.

=item B<-r> I<resampling_method>, B<--resampling=>I<resampling_method>

The resampling method to use. See gdalwarp(1) for the default and list
of possible values.

=item B<-b> I<width>[:I<color>], B<--border=>I<width>[:I<color>]

The I<width> and I<color> of the border to draw around each tile
(default: C<10:red>). I<width> is given in pixels, and I<color> needs to
comply to ImageMagick's Color Model Specification; see convert(1) for
details.
If one of I<width> or I<color> is left empty, the corresponding default
is used (but the other can be specified regardless, for instance in
C<:blue>).

=item B<-f> I<font>[:I<color>[:I<font>]], B<--font=>I<size>[:I<color>[:I<font>]]

The point I<size>, I<color> and I<font> to use when annotating each tile
with its name (default: C<25::Helvetica>). I<color> and I<font> need to
be ImageMagick-compatible; see convert(1) for details.
If I<color> is blank empty, the border's color is used instead.
If one of I<width> or I<font> is left empty, the corresponding default
is used (but the other can be specified regardless, for instance in
C<::DejaVu-Sans-Book>).

=item B<--s_srs=>I<SRS>

The source spatial reference set. If the input I<map>s include a SRS, it is
ignored. I<SRS> can be passed in any format format understood by
C<OGRSpatialReference.SetFromUserInput>, see gdalwarp(1) for details.

=item B<--t_srs=>I<SRS>

The target spatial reference set. This is required unless all I<map>s
have the same SRS, or if B<--s_srs> is used (in which case that SRS is
used for the target index).
I<SRS> can be passed in any format format understood by
C<OGRSpatialReference.SetFromUserInput>, see gdalwarp(1) for details.

=item B<--help>

Show a brief help.

=item B<--version>

Show the version numbers.

=item B<--man>

Show the manpage.

=back

=head1 EXAMPLES

=over 4

=item mkindex.pl *.tif miff:- | display -

Create an index with all TIFF files in the current directory; don't save
the index, but pipe it to display(1).

=item mkindex.pl -b '15:RGBA(255,0,0,.25)' --s_srs=epsg:3067 *.tif index.tif

Assume all map to be projected in EPSG:3067. Draw a 15-pixels wide , 25%
transparent red around each tile.

=item mkindex.pl -f '75::DejaVu-Sans-Book' --rename='s/U(.*)_RVK_25\.tif/$1/' *.tif index.tif

Use DejaVu Sans with a pointsize of 75 for the captions. Also, display
only the non-redundant part of the filenames.

=back

=head1 ENVIRONMENT

A temporary working directory is created under the B<TMPDIR> (or C</tmp>
if unset, at least in Unix:es), see C<File::Spec::tmpdir()>.

=head1 REQUIREMENTS

=over 4

=item *

Requires L<gdalwarp|http://gdal.org/> available on the command line.

=item *

Requires the perl modules
L<Inline|http://search.cpan.org/~sisyphus/Inline-0.53/C/C.pod>,
L<Geo::GDAL|http://geoinformatics.aalto.fi/doc/Geo-GDAL/html/>,
L<Geo::OSR|http://geoinformatics.aalto.fi/doc/Geo-GDAL/html/> and
L<Image::Magick|http://www.imagemagick.org/script/perl-magick.php>,
all available on L<CPAN|http://www.cpan.org/>.

=back

=head1 AUTHOR

Copyright 2013 Guilhem Moulin. mkindex is free software,
distributed under the GNU Public License, version 3 or later.

=cut

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 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   => sub {
                                # Geo::GDAL::TermProgress_nocb writes to STDOUT
                                open my $x, '>&', \*STDOUT or die "Can't dup STDOUT: $!";
                                open STDOUT, '>&', STDERR or die "Can't dup STDERR: $!";
                                Geo::GDAL::TermProgress_nocb(@_[0,1]);
                                open STDOUT, '>&', $x or die "Can't dup X: $!";
                                close $x;
                             }
             , debug      => 0
             , rename     => 's/\.[^.]+$//'
             , outsize    => [ 1024, 0 ]
             , border     => 10
             , bordercolor=> 'red'
             , fontsize   => 25
             , font       => 'Helvetica'
             , wm         => 512
);
my $gdal_opts = { INTERLEAVE => 'BAND'
                , TILED      => 'YES'
                , COMPRESS   => 'LZW'
                };

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 { pod2usage(2) unless $_[1] =~ /^(\d+)x(\d+)$/;
                                $config{outsize} = [ $1, $2 ];
                              }
          , "r|resampling=s"=> \$config{resampling}
          , "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);
# TODO: -n|--numeric: put numbers instead of text
# TODO: -o|--overview: use a different overview mapset
# TODO: when the input is RGB (not RGBA), we want add an alpha channel
# but keep the nodata value as matte within the map bounds.
# This seems to be possible with a clever use of
# gdalwarp/gdal_translate.
# http://www.gdal.org/frmt_gtiff.html
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;

    my $stripped;
    for (my $b=1; $b <= $ds->{RasterCount} && !$stripped; $b++) {
        my ($bx, $by) = $ds->GetRasterBand($b)->GetBlockSize();
        $stripped = 1 if ($bx == 1 and $by > 1) or ($bx > 1 and $by == 1);
    }
    warn "WARN: '$name' is stripped. Tiled structure is more efficient, especially for wrapping.\n"
        if $stripped;

    my $notband = $ds->GetMetadata('IMAGE_STRUCTURE')->{INTERLEAVE};
    undef $notband unless defined $notband and $notband ne 'BAND';
    warn "WARN: '$name' is '$notband'-interleaved; 'BAND' interleaving is recommended.\n"
        if $notband;

    if ($stripped or $notband) {
        warn "WARN: To fix '$name', run the following command in a POSIX-compatible shell.\n";
        warn "WARN: This is a lossless and usually size-preserving conversion.\n";
        my $old = $name;
        my ($new,$path,$suffix) = fileparse ($name, qr/\.[^.]*$/);
        $new = File::Spec->catfile($path,$new.'_new'.$suffix);
        my $driver = $ds->GetDriver()->{ShortName};
        my %opts = %$gdal_opts;
        my $compress = $ds->GetMetadata('IMAGE_STRUCTURE')->{COMPRESSION};
        $opts{COMPRESS} = $compress if $compress;
        warn "WARN:     ".join (' ',
             map { my $x = $_;
                   $x =~ s/"/\\"/;
                   $x = "\"$x\"" if $x =~ /[ ()';#{}*?~&|`!]/;
                   $x
                 }
                 ( 'gdal_translate', '-of', $driver
                 , (map { ('-co', $_.'='.$opts{$_}) } (keys %opts))
                 , $old, $new
                 )
             ). "\n";
    }

    my $srs;
    $map->{wkt} = $ds->GetProjection();
    unless ($config{s_srs}) {
        die "Error: Missing projection for '$name'.\n"
            unless $map->{wkt} or $config{s_srs};
        $map->{wkt} ||= $config{s_srs}->ExportToWkt();

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

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

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


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

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


# 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 = { BIGTIFF    => 'IF_NEEDED' };
my $driver_type = 'Byte';
die "Error: Unknown driver: '$driver'.\n"
    unless Geo::GDAL::GetDriverByName($driver);
$driver_opts->{$_} //= $gdal_opts->{$_} for (keys %$gdal_opts);

map &expand_rgba($_), @mapset;
my $mosaic = &merge();
&annotate( $mosaic );


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

# 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;
}

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

    &debug ("Opening '$map->{filename}'");
    my $src_ds = Geo::GDAL::Open($map->{filename}, 'ReadOnly') or exit 1;

    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);
    }

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

    $dst_ds->SetProjection( $map->{wkt} );
    $dst_ds->SetGeoTransform( $map->{s_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;
}


# Build the index. We use gdalwarp(1) because it does some clever
# computation to estimate the size of the reprojected tiles, and the
# GDALSuggestedWarpOutput() is not available in the Perl (or Python)
# API.
sub merge {
    my $dst_filename = &tempfile(SUFFIX => '.tif');
    &debug("Wrapping into file '$dst_filename'...");

    my @opts;
    # We use the optimisations suggested below
    # http://lists.osgeo.org/pipermail/gdal-dev/2009-August/021715.html
    # http://trac.osgeo.org/gdal/wiki/UserDocs/GdalWarp
    # '-multi' seems to slow down gdalwarp, though :-(
    push @opts, '--config', 'GDAL_CACHEMAX', $config{wm}, '-wm', $config{wm}
        if $config{wm};
    push @opts, '-q' unless $config{progress};
    push @opts, '-s_srs', $config{s_srs}->ExportToProj4() if $config{s_srs};
    push @opts, '-t_srs', $config{t_srs}->ExportToProj4() if $config{t_srs};
    push @opts, '-r', $config{resampling} if $config{resampling};
    push @opts, '-ts', $config{outsize}->[0], $config{outsize}->[1];
    my $driver_opts = $driver_opts;
    # gdalwarp cannot create properly compressed output
    # http://trac.osgeo.org/gdal/wiki/UserDocs/GdalWarp#GeoTIFFoutput-coCOMPRESSisbroken
    delete $driver_opts->{COMPRESS} if $driver_opts->{COMPRESS};
    push @opts, '-of', $driver, map { ('-co', $_.'='.$driver_opts->{$_}) }
                                    (keys %$driver_opts);
    my @src_filenames = map {$_->{filename}} @mapset;

    {
      # gdalwarp writes to STDOUT; we add >&2
      open my $x, '>&', \*STDOUT or die "Can't dup STDOUT: $!";
      open STDOUT, '>&', STDERR or die "Can't dup STDERR: $!";
      system 'gdalwarp', @opts, @src_filenames, $dst_filename;
      $? == 0 or die "gdalwarp failed.\n";
      open STDOUT, '>&', $x or die "Can't dup X: $!";
      close $x;
    }
    return $dst_filename;
}


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

    my $ds = Geo::GDAL::Open($filename, 'ReadOnly') or exit 1;
    my ($ok,@invt) = Geo::GDAL::InvGeoTransform([ $ds->GetGeoTransform() ]);
    undef $ds;
    die "Error: Couldn't inverse affine transformation\n" unless $ok;

    my $img = Image::Magick::->new();
    $img->Read( $filename );
    my ($width, $height) = $img->Get(qw/width height/);

    my @draw = ( fill => 'none'
               # list of supported colors: `convert -list color`
               , stroke => $config{bordercolor}
               , strokewidth => $config{border}
               , primitive => 'polygon'
               );
    my @ann = ( fill => $config{fontcolor} // $config{bordercolor}
              , strokewidth => 0
              # list of supported fonts: `convert -list font`
              , font => $config{font}
              , pointsize => $config{fontsize}
              );
    foreach my $map (@mapset) {
        # In case the projection have changed, the map may have been rotated,
        # so we draw a polygon instead of a rectangle. Also, we're using
        # a mask to avoid the text to be out of bounds.

        my @points = map {[ Geo::GDAL::ApplyGeoTransform(\@invt, @$_[0,1] ) ]}
                          @{$map->{t_corners}}[0,1,3,2];
        my $points = join (' ', map {join ',', @$_} @points);
        &debug( "Polygon: " . $points );
        # TODO: when drawing a semi-transparent border, all borders
        # should be in the same group so the corners are not "darker"
        # than the rest.
        $img->Draw( @draw, points => $points );

        # Find the middle of the polygon; We add a Y-offset to be
        # vertically centered.
        my ($x,$y) = (0,0);
        map { $x+= $_->[0] / (1.+$#points); $y+= $_->[1] / (1.+$#points); }
            @points;
        push @ann, x => $x, y => $y + $config{fontsize}/2., align=> 'Center';

        my $text = Image::Magick::->new( size => "${width}x${height}" );
        $text->ReadImage('xc:none');

        my $mask = Image::Magick::->new( size => "${width}x${height}" );
        $mask->ReadImage('xc:black');
        $mask->Draw( primitive => 'polygon'
                    , strokewidth => 0
                    , fill => 'white'
                    , points => $points );
        $text->Mask( mask => $mask );
        undef $mask;

        $text->Annotate( @ann, text => $map->{text} );
        $img->Composite( image => $text, compose => 'Over' );
        undef $text;
    }
    binmode STDOUT if $index =~ /^(?:[[:alnum:]]*:)?-$/; # Useful for pipes
    $img->Write($index);
    undef $img;
}


__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] ];
}