#!/usr/bin/env perl
# MMBLT to POV-Ray - Convert SwissTopo/L+T MMBL(T) DHM25 elevation data
# http://www.roe.ch/DHM25
# 
# Copyright (C) 2004-2005, Daniel Roethlisberger <daniel@roe.ch>
# 
# Redistribution and use, with or without modification, are permitted
# provided that the following conditions are met:
# 1. Redistributions must retain the above copyright notice, this list of
#    conditions and the following disclaimer.
# 2. The name of the author may not be used to endorse or promote products
#    derived from this software without specific prior written permission.
# 
# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
# OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
# NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
# THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
# 
# $Id: mmblt2pov.pl,v 1.3 2005/05/10 11:25:53 roe Exp $
use Getopt::Std;
use GD;
use strict;

use vars qw{ $opt_m $opt_p $opt_t $opt_o $opt_q $opt_h };
getopts("m:p:t:o:qh");
sub usage
{
    die "$0 [OPTIONS...] [INFILE]

Read a DHM25 matrix model in MMBLT/MMBL format from INFILE or stdin, and
convert it to a more universally useful PNG height map and POV-Ray object.

	-m file   Write PNG height map to file
	-p file   Write POV-Ray object source to file
	-t file   Use file as texture image for POV-Ray object
	-o XXXX   Equivalent of -m XXXX_hm.png -p XXXX.inc -t XXXX.png
	-q        Be quiet
	-h        This cruft

If none of -m, -p and -o are given, no output will be written, but the
MMBLT header and data will still be parsed and checked.

To transform mm1091.mlt into corresponding 1091_hm.png and 1091.inc, use:
	$0 -o 1091 mm1091.mlt
";
}
usage() if($opt_h);

my $heightmap_file = $opt_m || ($opt_o ? $opt_o .'_hm.png' : '');
my $povrayobj_file = $opt_p || ($opt_o ? $opt_o .'.inc' : '');
my $texture_file   = $opt_t || ($opt_o ? $opt_o .'.png' : '');

my $in_file = shift() || '-';

sub coord
{
	my $c = shift() || die;
	$c =~ s/\.0+$//;
	$c =~ s/([0-9]{2,3})([0-9]{3})/$1\/$2/;
	return $c;
}

if($in_file eq '-')
{
	*INFILE = \*STDIN;
}
else
{
	open(INFILE, "<$in_file") or die "Cannot open $in_file for reading!";
}

#NEWHEADER
#--------------------------------------------------------------------------------
#DHM25-MATRIXMODELL                             (c)BUNDESAMT F. LANDESTOPOGRAPHIE
#--------------------------------------------------------------------------------
#NORD-WEST ECKE     [M]  XXXX00.0  XXXX00.0     ERSTER HOEHENWERT
#SUED-OST ECKE      [M]  XXXX00.0  XXXX00.0     LETZTER HOEHENWERT
#MASCHENWEITE WE/NS [M]      25.0      25.0
#MATRIXDIMENSIONEN WE/NS    XXX       XXX       TOTAL    XXXXXX MATRIXPUNKTE
#HOEHENBEREICH     [DM]    XXXX      XXXX       (6 CHARACTER PRO HOEHENWERT)
#--------------------------------------------------------------------------------
#FORMAT                   ASCII                 L+T-FORMAT DHM25-MATRIXMODELL
#RECORDLAENGE(CHAR.)       XXXX                 XXX HOEHENWERTE PRO RECORD
#--------------------------------------------------------------------------------
#ENDHEADER
#  XXXX  XXXX  XXXX  [...]

my $header;
while(($_ = <INFILE>) && ($_ !~ /^NEWHEADER$/)) { }
die "No MMBLT header found!\n" unless($_);

$header = $_; # <- NEWHEADER
while(($_ = <INFILE>) && ($_ !~ /^ENDHEADER$/)) { $header .= $_; }
$header .= $_; # <- ENDHEADER

# FORMAT                   ASCII                 L+T-FORMAT DHM25-MATRIXMODELL
$header =~ /FORMAT +ASCII +L\+T-FORMAT DHM25-MATRIXMODELL/ or die "Not a DHM25 ASCII MMBLT file!";

# NORD-WEST ECKE     [M]  XXXX00.0  XXXX00.0     ERSTER HOEHENWERT
$header =~ /NORD-WEST ECKE +\[M\] +([0-9.]+) +([0-9.]+)/ or die "No N/W edge in MMBLT header!";
my $edge_x1 = $1;
my $edge_y1 = $2;

# SUED-OST ECKE      [M]  XXXX00.0  XXXX00.0     LETZTER HOEHENWERT
$header =~ /SUED-OST ECKE +\[M\] +([0-9.]+) +([0-9.]+)/ or die "No S/E edge in MMBLT header!";
my $edge_x2 = $1;
my $edge_y2 = $2;

# MASCHENWEITE WE/NS [M]      25.0      25.0
$header =~ /MASCHENWEITE WE\/NS \[M\] +([0-9.]+) +([0-9.]+)/ or die "No grid size in MMBLT header!";
my $grid_x = $1;
my $grid_y = $2;

# MATRIXDIMENSIONEN WE/NS    XXX       XXX       TOTAL    XXXXXX MATRIXPUNKTE
$header =~ /MATRIXDIMENSIONEN WE\/NS +([0-9]+) +([0-9]+) +TOTAL +([0-9]+) MATRIXPUNKTE/ or die "No dimensions in MMBLT header!";
my $dim_x = $1;
my $dim_y = $2;
my $total = $3;

# HOEHENBEREICH     [DM]    XXXX      XXXX       (6 CHARACTER PRO HOEHENWERT)
$header =~ /HOEHENBEREICH +\[DM\] +([0-9]+) +([0-9]+)/ or die "No height range in MMBLT header!";
my $h_min = $1;
my $h_max = $2;

die "Inconsistent header data ($dim_x * $dim_y != $total)!" if($dim_x * $dim_y != $total);

my $coord_x1 = coord($edge_x1);
my $coord_y1 = coord($edge_y1);
my $coord_x2 = coord($edge_x2);
my $coord_y2 = coord($edge_y2);

my $factor_x = ($dim_x - 1) * $grid_x / 1000;
my $factor_y = ($dim_y - 1) * $grid_y / 1000;
my $factor_h = ($h_max - $h_min) / 10000;
my $offset_x = $edge_x1 / 1000;
my $offset_y = $edge_y2 / 1000;
my $offset_h = $h_min / 10000;

print <<EOF unless($opt_q);
MMBLT Source File   : $in_file
Data Points in Model: $dim_x x $dim_y = $total
Data Grid Size   (m): $grid_x x $grid_y
Height Range    (dm): $h_min .. $h_max
Stretch Factors (km): $factor_x / $factor_y / $factor_h
N/W Edge Coordinates: $coord_x1//$coord_y1
S/E Edge Coordinates: $coord_x2//$coord_y2
EOF

my $img = new GD::Image($dim_x, $dim_y, 1) if($heightmap_file);
my $line;
for(my $y = 0; $y < $dim_y; $y++)
{
	for(my $x = 0; $x < $dim_x; $x++)
	{
		$line =~ s/^[\s\n\r]+//;
		$line = <INFILE> unless($line);
		die "Premature end of data!" unless($line);
		$line =~ /([0-9]+)/ or die "Illegal data!";
		$line = $';
		my $h = $1;
		die "Height value $h out of range $h_min..$h_max!" if($h < $h_min || $h > $h_max);
		my $h_n = ($h - $h_min) * 256 / ($h_max - $h_min);
		$img->setPixel($x,$y,$img->colorAllocate($h_n,$h_n,$h_n)) if($heightmap_file);
	}
}

if($in_file ne '-')
{
	close(INFILE);
}

if($heightmap_file)
{
	open(OUTFILE, ">$heightmap_file") or die "Cannot open $heightmap_file for writing!";
	print OUTFILE $img->png;
	close(OUTFILE);
}

if($povrayobj_file)
{
	open(OUTFILE, ">$povrayobj_file") or die "Cannot open $povrayobj_file for writing!";
	print OUTFILE <<EOF;
/*
 * POV-Ray Object Include File
 * Generated by $0
 * See http://www.roe.ch/DHM25
 * 
 * DHM25 Source File   : $in_file
 * Data Points in Model: $dim_x x $dim_y = $total
 * Data Grid Size   (m): $grid_x x $grid_y
 * Height Range    (dm): $h_min .. $h_max
 * Stretch Factors (km): $factor_x / $factor_y / $factor_h
 * N/W Edge Coordinates: $coord_x1//$coord_y1
 * S/E Edge Coordinates: $coord_x2//$coord_y2
 */

#declare FactorX = $factor_x;
#declare FactorY = $factor_y;
#declare FactorH = $factor_h;
#declare OffsetX = $offset_x;
#declare OffsetY = $offset_y;
#declare OffsetH = $offset_h;
#declare HeightMap = "$heightmap_file";
#declare TextureMap = "$texture_file";

object
{
 height_field
 {
  png HeightMap
  smooth
  rotate -90*x
 }
 pigment
 {
  image_map { png TextureMap }
 }
 scale <FactorX,FactorY,FactorH>
 translate <OffsetX,OffsetY,-OffsetH>
}
EOF
	close(OUTFILE);
}