2011-04-30 Map Drawing Using Polygons

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

use strict;
use CGI qw(:standard);
use SVG;
use Math::Geometry::Voronoi;
use Class::Struct;
use Math::Fractal::Noisemaker;
use List::Util qw(min max);
use Data::Dumper;

my $points = 3000;
my $width  = 1000;
my $height =  550;
my $center_x = $width / 2;
my $center_y = $height / 2;
my $radius = 500;

my %color = (beach => '#a09077',
	     ocean => '#44447a',);

struct World => { points => '@',
		  centroids => '@',
		  voronoi => '$',
		  height => '@',
		};

sub add_random_points {
  my ($world) = @_;
  for (my $i = 0; $i < $points; $i++) {
    push(@{$world->points}, [rand($width), rand($height)]);
  };
  # print(join("\n", map {join(",", $_->[0], $_->[1])} @{$world->points}));
  return $world;
}

sub add_voronoi {
  my ($world) = @_;
  $world->voronoi(Math::Geometry::Voronoi->new(points => $world->points));
  $world->voronoi->compute;
}

sub add_centroids {
  my ($world) = @_;
  $world->centroids([]); # clear
  foreach my $polygon ($world->voronoi->polygons) {
    push(@{$world->centroids}, centroid($polygon));
  }
}

sub centroid {
  my ($cx, $cy) = (0, 0);
  my $A = 0;
  my $polygon = shift;
  my ($point_index, @points) = @$polygon; # see Math::Geometry::Voronoi
  my $point = $points[$#points];
  my ($x0, $y0) = ($point->[0], $point->[1]);
  for $point (@points) {
    my ($x1, $y1) = ($point->[0], $point->[1]);
    $cx += ($x0 + $x1) * ($x0 * $y1 - $x1 * $y0);
    $cy += ($y0 + $y1) * ($x0 * $y1 - $x1 * $y0);
    $A += ($x0 * $y1 - $x1 * $y0);
    ($x0, $y0) = ($x1, $y1);
  }
  $A /= 2;
  $cx /= 6 * $A;
  $cy /= 6 * $A;
  return [$cx, $cy, $point_index];
}

sub add_height {
  my $world = shift;
  $Math::Fractal::Noisemaker::QUIET = 1;
  my $grid = Math::Fractal::Noisemaker::square();
  $world->height([]); # clear
  my $scale = max($height, $width); # grid is a square
  foreach my $point (@{$world->points}) {
    my $x = int($point->[0]*255/$scale);
    my $y = int($point->[1]*255/$scale);
    my $h = 0; # we must not skip any points!
    $h = $grid->[$x]->get($y) / 255
      unless $x < 0 or $y < 0 or $x > 255 or $y > 255;
    push(@{$world->height}, $h);
  }
}

sub raise_point {
  my ($world, $x, $y, $radius) = @_;
  my $i = 0;
  foreach my $point (@{$world->points}) {
    my $dx = $point->[0] - $x;
    my $dy = $point->[1] - $y;
    my $d = sqrt($dx * $dx + $dy * $dy);
    my $v = max(0, $world->height->[$i] - $d / $radius);
    $world->height($i, $v);
    $i++;
  }
}

sub svg {
  my $world = shift;
  my $svg = new SVG(-width => $width,
		    -height => $height, );
  foreach my $polygon ($world->voronoi->polygons) {
    my ($point_index, @points) = @$polygon; # see Math::Geometry::Voronoi
    my $x = $world->points->[$point_index]->[0];
    my $y = $world->points->[$point_index]->[1];
    next if $x < 0 or $y < 0 or $x > $width or $y > $height;
    my $z = int($world->height->[$point_index] * 255);
    my $color = $z == 0 ? $color{ocean} : "rgb($z,$z,$z)";
    my $path = join(",", map { map { int } @$_ } @points);
    $svg->polygon(points => $path,
		  fill => $color,
		  style => { 'stroke-width' => 1,
			     'stroke' => 'black'});
  }
  return $svg->xmlify();
}

sub response {
  print header(-type=>'image/svg+xml');
  print shift;
}

sub main {
  if (path_info eq '/source') {
    seek DATA, 0, 0;
    print "Content-type: text/plain; charset=UTF-8\r\n\r\n", <DATA>;
  } else {
    srand(param('seed') || time);
    my $world = new World;
    add_random_points($world, $points);
    add_voronoi($world);
    for (my $i = 2; $i--; ) {
      # Lloyd Relaxation
      add_centroids($world);
      $world->points($world->centroids);
      add_voronoi($world);
    }
    # skip corner improvement
    # skip Delaunay triangulation
    add_height($world);
    raise_point($world, $center_x, $center_y, $radius);
    # draw
    response(svg($world));
  }
}

main ();

__DATA__