#!/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 warnings;
use strict;

=head1 NAME

mkindex.pl - Create an index file from a mapset

=head1 SYNOPSIS

B<mkindex.pl> [B<-q>] [B<-d>] [B<--rename> I<perlexpr>] [B<-n> I<file>]
[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>] [-o I<map>] [B<--crop>] 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<-n> I<file>, B<--numeric>=I<file>

Instead of annoting the index with the map names, use an increasing
counter starting from 1, and write the mapping I<num: text> in I<file>,
where I<text> is a map name massaged with I<perlexpr> (see B<--rename>).

=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<-o> I<map>, B<--overview=>I<map>

Use I<map> as the overview map. This option can be used serveral times
to use multiple (warped) overview maps.

=item B<--crop>

Alternative overview map may be much bigger than the actual mapset. Use
B<--crop> to trim away the useless part. Note that the resulting index file
may be smaller than wanted.

=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 '25:RGBA(255,0,0,.5)' -f ':RGBA(0,0,255,.66)' --s_srs=epsg:3067 *.tif index.tif

Assume all map to be projected in EPSG:3067. Draw a 25-pixels wide, 50%
transparent red border around each tile, and write the caption in 66%
opaque blue.

=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.

=item mkindex.pl -o overview1.tif -o overview2.tif --crop maps*.tif index.tif

Use the C<overview1.tif> and C<overview2.tif> files as overview maps for
the mapset C<maps*.tif>. Also, crop the parts of the index that are not
in the mapset.

=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>.
They are 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 File::Spec;
use POSIX qw/floor ceil/;
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 "Error: Can't dup STDOUT: $!";
                                open *STDOUT, '>&', \*STDERR or die "Error: Can't dup STDERR: $!";
                                Geo::GDAL::TermProgress_nocb(@_[0,1]);
                                open *STDOUT, '>&', $x or die "Error: 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]
                              }
          , "n|numeric=s" => \$config{numeric}
          , "s|size=s" => sub { pod2usage(2) unless $_[1] =~ /^(\d+)x(\d+)$/;
                                $config{outsize} = [ $1, $2 ];
                              }
          , "r|resampling=s"=> \$config{resampling}
          , "o|overview=s@" => \$config{overviews}
          , "crop"          => \$config{crop}
          , "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 { my %h = ( WKT => Geo::OSR::GetUserInputAsWKT($_[1]) );
                               $config{s_srs} = Geo::OSR::SpatialReference::->new(%h); }
          , "t_srs=s" => sub { my %h = ( WKT => Geo::OSR::GetUserInputAsWKT($_[1]) );
                               $config{t_srs} = Geo::OSR::SpatialReference::->new(%h); }
          , "man"     => sub { pod2usage(-exitstatus => 0, -verbose => 2) }
          )
    or pod2usage(2);
# 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;
$config{fontcolor} //= $config{bordercolor};
if ($config{numeric}) {
    $config{numeric_format} = "%.".(length $#mapset)."d";
    open my $fh, '>', $config{numeric}
        or die "Error: Can't open '".$config{numeric}."': $!";
    $config{numeric} = $fh;
}

&debug( "Using mapset:", map {"  - ".$_->{filename}} @mapset );
if ($config{overviews}) {
    $config{overviews} = [ map {{ filename => $_ }} @{$config{overviews}} ];
    &debug( "Using overview:", map {"  - ".$_->{filename}} @{$config{overviews}} );
}
else {
    $config{overviews} = \@mapset;
}
&debug( "Building index: $index" );


# Find the source and target SRS, and ensure the source files are in a
# decent format.
my $s_srs;
foreach my $map (@{$config{overviews}}) {
    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', paste 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,$ext) = fileparse ($name, qr/\.[^.]*$/);
        $new = File::Spec->catfile($path,$new.'_new'.$ext);
        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";
    }

    if ($config{s_srs}) {
        # The 's_srs' option takes precedence.
        $map->{wkt} = $config{s_srs}->ExportToWkt();
    }
    else {
        # Croak if the source files don't have a valid projection.
        $map->{wkt} = $ds->GetProjection()
            or die "Error: Missing projection for '$name'.\n";

        # Croak unless either 't_srs' is set, or all the map are in the
        # same projection.
        my $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;
        }
    }
}

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


# Store the corners of each tile, an the text to display.
my $num = 0;
foreach my $map (@mapset) {
    my $name = $map->{filename};

    &debug ("Opening '$name'");
    my $ds = Geo::GDAL::Open($name, 'ReadOnly') or exit 1;

    my $geot = [ $ds->GetGeoTransform() ];
    $map->{t_corners} = [ map [ Geo::GDAL::ApplyGeoTransform($geot, @$_) ],
                              ( [ 0,                    0 ]
                              , [ 0,                    $ds->{RasterYSize}-1 ]
                              , [ $ds->{RasterXSize}-1, 0 ]
                              , [ $ds->{RasterXSize}-1, $ds->{RasterYSize}-1 ]
                              )
                          ];

    my $srs = $config{s_srs};
    unless ($srs) {
        my $wkt = $ds->GetProjection()
            or die "Error: Missing projection for '$name'.\n";
        $srs = Geo::OSR::SpatialReference::->new( $wkt );
    }
    &debug( "Using source SRS: ". $srs->ExportToProj4() );

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

        $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} = $_; }
    if ($config{numeric}) {
        my $fh = $config{numeric};
        printf $fh $config{numeric_format}.': '.$map->{text}."\n", ++$num;
        $map->{text} = sprintf $config{numeric_format}, $num;
    }
}
close $config{numeric} if $config{numeric};


# 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);
my $driver_ext = Geo::GDAL::GetDriverByName($driver)->Extension();

&expand_rgba($_) foreach @{$config{overviews}};
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;
    my $name = $map->{filename};

    &debug ("Opening '$name'");
    my $src_ds = Geo::GDAL::Open($name, '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.");
    $map->{filename} = &tempfile(SUFFIX => '.'.$driver_ext);
    if ($config{progress}) {
        print STDERR "Processing input file $name.\n";
        &{$config{progress}}(0.0);
    }

    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 ($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([ $src_ds->GetGeoTransform() ]);
    $dst_ds->SetGCPs( $src_ds->GetGCPs(), $src_ds->GetGCPProjection() )
        if $src_ds->GetGCPCount();

    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++) {
            # 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/";
            my $dst_data = &subst( $src_data, $lookup->[$b], $xSize );
            $dst_ds->Band(1+$b)->WriteRaster(0, $y, $xSize, 1, $dst_data);
        }

        &{$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 => '.'.$driver_ext);
    &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}, @{$config{overviews}};
    {
      # gdalwarp writes to STDOUT; we add >&2
      open my $x, '>&', \*STDOUT or die "Error: Can't dup STDOUT: $!";
      open *STDOUT, '>&', \*STDERR or die "Error: Can't dup STDERR: $!";
      system 'gdalwarp', @opts, @src_filenames, $dst_filename;
      $? == 0 or die "Error: gdalwarp failed.\n";
      open *STDOUT, '>&', $x or die "Error: Can't dup X: $!";
      close $x;
    }
    return $dst_filename;
}


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

    # Inverse the transformation, so that we can convert from projected
    # coordinates to pixel coordinates.
    my $ds = Geo::GDAL::Open($filename, 'ReadOnly') or exit 1;
    my @invt = Geo::GDAL::InvGeoTransform([ $ds->GetGeoTransform() ]);
    undef $ds;

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

    my @draw = ( fill => 'none'
               , strokewidth => $config{border}
               , primitive => 'polygon'
               );
    my @ann = ( stroke => 'none'
              # To get the list of supported fonts: `convert -list font`.
              , font => $config{font}
              , pointsize => $config{fontsize}
              , fill => $config{fontcolor}
              , gravity => 'Center'
              );

    # Separate RGB and alpha channels. This complicated business is
    # because we want all borders to share the same alpha channel
    # instead of having more opaque intersections of semi-transparent
    # areas (due to the default composition, which multiplies the
    # values.
    #
    # To get the list of known color names: `convert -list color`.
    my $pix = Image::Magick::->new();
    my ($r, $g, $b, $a) = $pix->QueryColor( $config{bordercolor} );
    my $max = (2 << ($pix->Get('depth') - 1)) - 1.;
    undef $pix;
    ($r,$g,$b,$a) = map { ($_ // 0) * 100. / $max } ($r,$g,$b,$a);
    my ($b_rgb, $b_alpha) = ("RGB($r%,$g%,$b%)", 'GRAY('.(100-$a).'%)');

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

    my $mask;
    if ($a > 0) {
        &debug( "Adding a transparency mask (border opacity: $a%)." );
        $mask = Image::Magick::->new( size => "${width}x${height}" );
        $mask->Read( 'xc:black' );
        # We really want 'Gray', but a non-linear Gray. See
        # http://www.imagemagick.org/script/color-management.php
        # for the trick:
        #   convert img.png -set colorspace RGB -colorspace gray RGBimg.png
        # (Also, setting the 'colorspace' during the object creation was
        # useless.)
        $mask->Set( colorspace => 'RGB' );
    }

    my ($lx0,$uy0,$rx0,$ly0);
    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.

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

        # Add the borders. We do that after the text to draw over in
        # case of intersection.
        $rgb->Draw( @draw, stroke => $b_rgb,   points => $points );
        $mask->Draw(@draw, stroke => $b_alpha, points => $points ) if $mask;

        # Get the extremes of this polygon.
        my ($lx,$uy,$rx,$ly);
        foreach (@points) {
            $lx = floor $_->[0] unless defined $lx and $lx < floor $_->[0];
            $uy = floor $_->[1] unless defined $uy and $uy < floor $_->[1];
            $rx = ceil  $_->[0] unless defined $rx and $rx > ceil  $_->[0];
            $ly = ceil  $_->[1] unless defined $ly and $ly > ceil  $_->[1];
        };
        # The region where to confine the caption.
        my $cell = ($rx-$lx).'x'.($ly-$uy).'+'.$lx.'+'.$uy;
        &debug( "Cell: $cell" );

        # Add the caption. ('gravity' is relative to the 'geometry'.)
        # Everything that's not in the cell's region is trimmed away.
        $raster->MogrifyRegion( $cell, 'Annotate'
                              , @ann
                              , geometry => $cell
                              , text => $map->{text} );

        if ($config{crop}) {
            $lx0 = $lx unless defined $lx0 and $lx0 < $lx;
            $uy0 = $uy unless defined $uy0 and $uy0 < $uy;
            $rx0 = $rx unless defined $rx0 and $rx0 > $rx;
            $ly0 = $ly unless defined $ly0 and $ly0 > $ly;
        }
    }

    # Change the colorspace (now to non-linear gray), and turn off
    # the alpha channel. Then copy this (shared) channel back into the
    # RGB annotation channel.
    if ($mask) {
        $mask->Set( colorspace => 'Gray', alpha => 'Off' );
        $rgb->Composite( image => $mask, compose => 'CopyOpacity' );
        undef $mask;
    }

    # Merge all the annotations into the mosaic.
    $raster->Composite( image => $rgb, compose => 'Over' );
    undef $rgb;

    # Throw away the outside.
    if ($config{crop}) {
        my $crop = ($rx0-$lx0).'x'.($ly0-$uy0).'+'.$lx0.'+'.$uy0;
        &debug( "Crop: $crop" );
        $raster->Crop( geometry => $crop );
    }

    binmode STDOUT if $index =~ /^(?:[[:alnum:]]*:)?-$/; # Useful for pipes
    $raster->Write($index);
    undef $raster;
}

__END__

__C__

SV *subst(unsigned char *data, unsigned char *lookup, int length) {
    char *out = malloc(length * sizeof(unsigned char));
    for (int i = 0; i < length; i++)
        out[i] = lookup[ data[i] ];
    SV* ret = newSVpvn(out, length);
    free(out);
    return(ret);
}