#!/usr/bin/perl


BEGIN {
  push @INC , '../';  
  push @INC , '../SVG';
}

use SVG;

#use Getopt::Long;

my $text = 0;
my $smooth = 2;
my $max = 50;
my $min = 20;
my $grid = $min+int(rand($max-$min));
my $scale=200/$grid;

#GetOptions(
#           "t|text" => \$text,
#           "s|smooth=i" => \$smooth,
#           "r|res=i" => \$grid,
#);



my @points;
my @fixed;
my $temp;
my $index = 0;
my $lowest=1000;
my $highest=-1000;

my @colour_map;
my $image;
my $nodes;

#
# read the data and build an array for the fixed points
#


my @data = (
"'".int($grid/2).','.int($grid/2).",'",
#"'".int(rand($max-$min)).','.int(rand($max-$min)).','.int(rand(100))."'",
"'".int(0.7*$grid).','.int(0.7*$grid).",100'",
);
foreach my $i (@data) {
#while (<DATA>) {
  chomp $i;
  if ($i) {
    my ($x,$y,$temp) = split (',',$i);
    $points[$x][$y] = $temp;
    $fixed[$index][0] = $x;
    $fixed[$index][1] = $y;
    $fixed[$index][2] = $temp;
    $nodes = $index;
    $index++;
  }
}

#
#Then I loop over the whole grid and interpolate each point that has
#temperature zero (ie. uninitialized by a measurement)... 
#

#
# interpolate cells with no temperature
#

for my $x (0..$grid) {
  for my $y (0..$grid) {
    if (undef $points[$x][$y]) {$points[$x][$y] = 0};
    #if ($points[$x][$y] == 0) {
      $points[$x][$y] = &interpolate($x,$y,\@fixed, $smooth,$nodes);
    #}

    if ($points[$x][$y] > $highest) { $highest = $points[$x][$y]; }
    if ($points[$x][$y] < $lowest) { $lowest = $points[$x][$y]; }
    print "temp ($x,$y) = $points[$x][$y]\n" if ($debug);
  }
}

#
#Finally, if the text flag wasn't thrown, I build a graphics image using the GD module, 
#putting the output on stdout for ease of use in CGI if necessary: 
#

if (! $text) {

    #
    # initialize the GD image object
    #

    $image = SVG->new(width=>$scale*($grid+50),height=>$scale*($grid+50));

    #
    # build a nice colour map
    #

    @colour_map = &mkmap;

    # 
    # draw the points, temps as colours
    #
    $image->title(id=>'title')->cdata('A random surface map in SVG');
    $image->desc(id=>'title')->cdata('This SVG is being generated using SVG.pm V 1.12 and a set of random numbers. Sorry it takes so long to run, I am not tarring the data.');
    $image->text(x=>80,y=>20)->cdata("A $grid x $grid thermal map");
    $image->comment(-comment=>'Draw the image map.');  
    for my $x (0..$grid) {
        for my $y (0..$grid) {
            
            my $index = int( (($points[$x][$y]-$lowest)/($highest-$lowest+0.0001)) * 199 );
            $image->rect(x=>$x*$scale+25,y=>$y*$scale+25,width=>$scale,
                    height=>$scale,
                    style=>{fill=>"$colour_map[$index]",stroke=>"$colour_map[$index]"
                    });

            print "index - $index -- $colour_map[$index]\n" if ($debug);
        }
    }

    #
    # draw a strip with the colour map in it
    #
    $image->comment('Draw the strip');
    for my $i (0..($grid+1)) {
        $image->rect(x=>($grid+6)*$scale+25,y=>$i*$scale+25,                 
                width=>5*$scale,height=>$scale,
                style=>{fill=>$colour_map[199-($i/($grid+1))*199]});
    }

    #
    # and output as SVG
    #
    print "Content-type: image/svg+xml\n\n";
    
    print $image->render;
}


# Of course, the important bits are in the subroutines. interpolate sets up... 
# subroutines -------------------------------------------------

sub interpolate {
  my $x = shift;
  my $y = shift;
  my $fixed = shift;
  my $smoothness = shift;
  my $nodes = shift;

  if (! defined($smoothness)) { $smoothness = 2; }

  my @distances;
  my $coefs;
  my $cosum;
  my $i;
  my $temp;

#...and then does the real work. For the point to interpolate, we look at each fixed measurement and calculate the distance. Note that if your measured points never change coordinates,
#you could pre-calculate a table of coefficients and just look them up, which would be a zillion times faster. 

  for $i (0..$nodes) {
    

    $distance[$i] = sqrt( ($x - $fixed->[$i]->[0])**2 +
                          ($y - $fixed->[$i]->[1])**2
                        );
    print "x='$x' dx='$fixed->[$i]->[0]', y= '$y',dy='$fixed->[$i]->[1]' distance =  $distance[$i]\n" if ($debug);

  }

#Then we build some co-efficients that are inverse squares of the distance, and add up the known temperatures times their distance coefficients... 

  for $i (0..$nodes) {
    $coefs = eval(1/(0.001+$distance[$i]**$smoothness));
    $cosum += $coefs;
    $temp += ($coefs * $$fixed[$i][2]);
  }

#Notice we saved the sum of the coefficients? That's because now we divide by the sum of coefficients to get a temperature that is an average of the fixed temperatures, where the fixed
#temperatures are weighted by distance. Clever! Notice we also have a little exception to dodge divide by zero events. 

  unless ($cosum) { $cosum = 1 };
  $temp /= $cosum;

  return $temp;
  print "temp = $temp\n" if ($debug);
}

#And the mkmap routine just builds a funky colour map. 


sub mkmap {

  # build a nice colour map for the display.  Note that in this map
  # pure green (0,255,0) has two entries.  This has the effect of
  # highlighting the contour at the mean with a bright yellow line.
  # I admit, it was a fencepost error that turned out to be a
  # feature.

  my @colour_map;
  $colour_map[255] = $image->colorAllocate(0,0,0);

  for my $i (0..99) {
    my $red = int(($i/100) * 255);
    my $green = int(($i/100) * 255);
    my $blue = int((1-($i/100)) * 255); 
    $colour_map[$i] = $image->colorAllocate($red,$green,$blue);
  }

  for (0..99) {
    my $red = 255;
    my $green = int((1-($_/100)) * 255); 
    my $blue = 0;

    $colour_map[$_+100] = $image->colorAllocate($red,$green,$blue);
  }
    
  return @colour_map;
}