#!/usr/local/bin/perl 
use utf8;
use 5.010;

use Data::Dumper;
use GD;
use Math::Random;
use Getopt::Long;
use Pod::Usage;
use File::Basename;

my $man = 0;
my $help = 0;

#
# Default values for options:
#
my $out_width = 3600;
my $out_height = 2550;
my $patch_dir;

#
# Parse options.
#
my $getopt_result = GetOptions (
  "help|?" => \$help, "man" => \$man,
  "height=i" => \$out_height,
  "width=i" => \$out_width,
  "patch_dir=s" => \$patch_dir,
  "verbose+" => \$verbose) or pod2usage(2);
pod2usage(1) if $help;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;

my $out_name = shift @ARGV;
if (!defined($out_name) or $out_name eq "") {
  die "output_file is required; use --help for help.\n";
}
my $out_image = GD::Image->new($out_width, $out_height,1);

#
# Read in each patch file and place it in the patches hash.
# The keys to the patches hash are the heights of the patches.
# The height of each patch is its file name.
#
my %patches;
my @patch_files = glob "${patch_dir}/*.png";
foreach my $patch_file (@patch_files) {
  if ($verbose > 1) {
    print "Opening ", $patch_file, ".\n";
  }
  if (! (-r $patch_file)) {
    die "File ", $patch_file, " cannot be read.\n";
  }
  my $image_data = GD::Image->newFromPng($patch_file, 1);
  my ($filename, $directory, $extension) = fileparse ($patch_file, ".png");
  my $patch = ($filename + 0);
  $patches{$patch} = [ $patch_file, $image_data ];
}

#
# Compute the mean and standard deviation of each color
# component for each patch.
# We use the scRGB color system, which is a 13-bit linear coding
# based on sRGB's 8-bit companded coding.
#
foreach my $patch (keys %patches) {
  my ($image_name, $image_data)  = @{$patches{$patch}};
  if ($verbose > 1) {
    print $patch, " => ", $image_name, "\n";
  }
  my $image_width = $image_data->width;
  my $image_height = $image_data->height;
  my $pixel_count = $image_width * $image_height;
  my $red_sum = 0;
  my $green_sum = 0;
  my $blue_sum = 0;
  for (my $y = 0; $y < $image_height; $y++) {
    for (my $x = 0; $x < $image_width; $x++) {
      my $index = $image_data->getPixel($x, $y);
      my ($red_val, $green_val, $blue_val) = $image_data->rgb($index);
      $red_val = &to_scrgb ($red_val);
      $green_val = &to_scrgb ($green_val);
      $blue_val = &to_scrgb ($blue_val);
      $red_sum += $red_val;
      $green_sum += $green_val;
      $blue_sum += $blue_val;
    }
  }
  my $red_mean = $red_sum / $pixel_count;
  my $green_mean = $green_sum / $pixel_count;
  my $blue_mean = $blue_sum / $pixel_count;

  my $red_dev = 0;
  my $green_dev = 0;
  my $blue_dev = 0;
  for (my $y = 0; $y < $image_height; $y++) {
    for (my $x = 0; $x < $image_width; $x++) {
      my $index = $image_data->getPixel($x, $y);
      my ($red_val, $green_val, $blue_val) = $image_data->rgb($index);
      $red_val = &to_scrgb ($red_val);
      $green_val = &to_scrgb ($green_val);
      $blue_val = &to_scrgb ($blue_val);
      my $red_deviation = ($red_mean - $red_val);
      $red_dev += ($red_deviation)**2;
      my $green_deviation = ($green_mean - $green_val);
      $green_dev += ($green_deviation)**2;
      my $blue_deviation = ($blue_mean - $blue_val);
      $blue_dev += ($blue_deviation)**2;
    }
  }
  my $red_stddev = sqrt($red_dev / $pixel_count);
  my $green_stddev = sqrt($green_dev / $pixel_count);
  my $blue_stddev = sqrt($blue_dev / $pixel_count);

#
# Augment the hash value with the statistics.
#
  $patches{$patch} = [ $image_name, $image_data, $red_mean, $green_mean, $blue_mean,
    $red_stddev, $green_stddev, $blue_stddev ];
  if ($verbose > 1) {
    print "  (", $red_mean,  "[", $red_stddev, "], ", 
      $green_mean, "[", $green_stddev, "], ",
      $blue_mean, "[", $blue_stddev, "]).\n";
  }
}

#
# If --patch_dir was not specified, construct a default sky.
#
my @patch_heights = sort {$a <=> $b} keys %patches;
my $patch_count = @patch_heights;
if ($patch_count == 0) {
  my %model_patches = (
     100 => ["Default_01", undef, 409, 1447, 3037, 110, 93, 108],
     300 => ["Default_02", undef, 458, 1579, 3212, 120, 98, 105],
     500 => ["Default_03", undef, 535, 1725, 3411, 130, 104, 109],
     700 => ["Default_04", undef, 644, 1920, 3643, 136, 108, 115],
     900 => ["Default_05", undef, 761, 2149, 3901, 139, 115, 124],
    1100 => ["Deafult_06", undef, 913, 2418, 4207, 148, 126, 128],
    1300 => ["Default_07", undef, 1135, 2739, 4502, 164, 142, 121],
    1500 => ["Default_08", undef, 1456, 3132, 4800, 186, 153, 122],
    1700 => ["Deafult_09", undef, 1898, 3540, 5065, 197, 151, 108],
    1900 => ["Deafult_10", undef, 2417, 3930, 5224, 212, 133, 91],
    2000 => ["Default_11", undef, 2678, 4071, 5249, 197, 114, 88],
    2100 => ["Default_12", undef, 2841, 4121, 5212, 147, 96, 109],
    2212 => ["Default_13", undef, 2800, 3972, 4948, 170, 175, 174],
  );
#
# The model patches assume a total height of 2319.  Scale the defaults based
# on the actual height.
#
  foreach my $patch (sort {$a <=> $b} keys %model_patches) {
    my $new_height = &roundoff ($patch * ($out_height / 2319));
    if ($verbose > 1) {
      print "Height ", $patch, " becomes ", $new_height, ".\n";
    }
    @patches{$new_height} = @model_patches{$patch};
  }
  @patch_heights = sort {$a <=> $b} keys %patches;
  $patch_count = @patch_heights;
}

if ($verbose > 0) {
  print "Patch statistics: \n";
  foreach my $patch (sort {$a <=> $b} keys %patches) {
    my ($image_name, $image_data, $red_mean, $green_mean, $blue_mean,
      $red_stddev, $green_stddev, $blue_stddev) = @{@patches{$patch}};
    print "    ", $patch, " => [\"", $image_name, "\", undef, ", &roundoff ($red_mean),
      ", ", &roundoff ($green_mean), ", ", &roundoff ($blue_mean), ", ",
    &roundoff ($red_stddev), ", ", &roundoff ($green_stddev), ", ", 
    &roundoff ($blue_stddev), "],\n";
  }
}

#
# Interpolate and extrapolate from the patches we are given to every
# row of the output image.  If there are no patches, we construct the
# default sky.  If there is only one patch, all rows
# get its statistics.  If there are two, we construct a line through
# the two points for each statistic, and all rows are on that line.
# If there are more than two we do a piecewise-linear line.  The first
# and last lines extend to the top and bottom of the image.
#
my @patch_heights = sort {$a <=> $b} keys %patches;
my $patch_count = @patch_heights;

if ($patch_count == 1) {
  my $patch = $patch_heights [0];
  my ($image_name, $image_data, $red_mean, $green_mean, $blue_mean,
    $red_stddev, $green_stddev, $blue_stddev) = @{@patches{$patch}};
  if ($verbose > 0) {
    print "Only segment: (", $red_mean, ", ", $green_mean, ", ", $blue_mean, ").\n";
  }
  for (my $y = 0; $y < $out_height; $y++) {
    @patches{$y} = [ $image_name, $image_data, $red_mean, $green_mean,
      $blue_mean, $red_stddev, $green_stddev, $blue_stddev ] unless defined ($patches{$y});
  }
}

if ($verbose > 2) {
  print Dumper(%patches);
}

for (my $patch_index = 0; $patch_index < ($patch_count - 1); $patch_index++) {
  my $patch_high = $patch_heights [$patch_index];
  my $patch_low = $patch_heights [$patch_index + 1];
  my $high_height;
  my $low_height;
  if ($patch_index > 0) {
    $high_height = $patch_high;
  } else {
    $high_height = 0;
  }
  if ($patch_index < ($patch_count - 2)) {
    $low_height = $patch_low;
  } else {
    $low_height = $out_height;
  }
  if ($verbose > 1) {
    print "Segment ", $patch_index, " from ", $patch_high, " to ", $patch_low, "\n";
  }
  my ($image_name_high, $image_data_high, $red_mean_high, $green_mean_high,
    $blue_mean_high, $red_stddev_high, $green_stddev_high, $blue_stddev_high) =
    @{@patches{$patch_high}};
  my ($image_name_low, $image_data_low, $red_mean_low, $green_mean_low,
    $blue_mean_low, $red_stddev_low, $green_stddev_low, $blue_stddev_low) =
    @{@patches{$patch_low}};
  if ($verbose > 1) {
    print "  (", $red_mean_high, ", ", $green_mean_high, ", ", $blue_mean_high, ")\n";
    print "  (", $red_mean_low, ", ", $green_mean_low, ", ", $blue_mean_low, ")\n";
  }
  my $red_mean_slope = ($red_mean_low - $red_mean_high) / ($patch_low - $patch_high);
  my $red_mean_intercept = $red_mean_high - ($patch_high * $red_mean_slope);
  my $red_stddev_slope = ($red_stddev_low - $red_stddev_high) / ($patch_low - $patch_high);
  my $red_stddev_intercept = $red_stddev_high - ($patch_high * $red_stddev_slope);
  my $green_mean_slope = ($green_mean_low - $green_mean_high) / ($patch_low - $patch_high);
  my $green_mean_intercept = $green_mean_high - ($patch_high * $green_mean_slope);
  my $green_stddev_slope = ($green_stddev_low - $green_stddev_high) / ($patch_low - $patch_high);
  my $green_stddev_intercept = $green_stddev_high - ($patch_high * $green_stddev_slope);
  my $blue_mean_slope = ($blue_mean_low - $blue_mean_high) / ($patch_low - $patch_high);
  my $blue_mean_intercept = $blue_mean_high - ($patch_high * $blue_mean_slope);
  my $blue_stddev_slope = ($blue_stddev_low - $blue_stddev_high) / ($patch_low - $patch_high);
  my $blue_stddev_intercept = $blue_stddev_high - ($patch_high * $blue_stddev_slope);

  if ($verbose > 1) {
    print $high_height, ":", $low_height, " -> (", $red_mean_slope, ", ", $green_mean_slope,
      ", ", $blue_mean_slope, ")\n";
    print "  (", $red_mean_intercept, ", ", $green_mean_intercept, ", ", $blue_mean_intercept, 
      ")\n";
  }
  for (my $y = $high_height; $y < $low_height; $y++) {
    my $red_mean = ($y * $red_mean_slope) + $red_mean_intercept;
    my $red_stddev = ($y * $red_stddev_slope) + $red_stddev_intercept;
    my $green_mean = ($y * $green_mean_slope) + $green_mean_intercept;
    my $green_stddev = ($y * $green_stddev_slope) + $green_stddev_intercept;
    my $blue_mean = ($y * $blue_mean_slope) + $blue_mean_intercept;
    my $blue_stddev = ($y * $blue_stddev_slope) + $blue_stddev_intercept;
    if ($verbose > 1) {
      printf "y = %g, rgb = (%7.3f, %7.3f, %7.3f), std. dev. = (%7.3f, %7.3f, %7.3f)\n", $y,
        $red_mean, $green_mean, $blue_mean, $red_stddev, $green_stddev, $blue_stddev;
    }
    @patches{$y} = [ $image_name_high, $image_data_high, $red_mean, $green_mean,
      $blue_mean, $red_stddev, $green_stddev, $blue_stddev ] unless defined ($patches{$y});
  }
}

#
# Now generate the image.
#
if ($verbose > 0) {
  print "Generating a ", $out_width, " by ", $out_height, " image.\n";
}
for (my $y = 0; $y < $out_height; $y++) {
 my ($image_name, $image_data, $red_mean, $green_mean, $blue_mean,
      $red_stddev, $green_stddev, $blue_stddev) = @{@patches{$y}};

  if ($verbose > 1) {
    printf "y = %g, rgb = (%7.3f, %7.3f, %7.3f), std. dev. = (%7.3f, %7.3f, %7.3f)\n", $y, 
      $red_mean, $green_mean, $blue_mean, $red_stddev, $green_stddev, $blue_stddev;
  }
  for (my $x = 0; $x < $out_width; $x++) {
    my $red_val = &random_component_value ($red_mean, $red_stddev);
    my $green_val = &random_component_value ($green_mean, $green_stddev);
    my $blue_val = &random_component_value ($blue_mean, $blue_stddev);
    $red_val = &to_srgb ($red_val);
    $green_val = &to_srgb ($green_val);
    $blue_val = &to_srgb ($blue_val);
    $red_val = &roundoff ($red_val);
    $green_val = &roundoff ($green_val);
    $blue_val = &roundoff ($blue_val);
    my $pixel_color = $out_image->colorAllocate ($red_val, $green_val, $blue_val);
    $out_image->setPixel($x, $y, $pixel_color);
  }
}

open OUTFILE, ">", $out_name;
binmode (OUTFILE);
print OUTFILE $out_image->png;
close OUTFILE;

#
# Given a mean and standard deviation, generate a random
# value for a color component.  The values are clamped
# to (0,8192) to keep them within the sRGB color gamut.
#
sub random_component_value {
  my ($mean, $std_dev) = @_;
  my $value = random_normal (1, $mean, $std_dev);
  $value = 0 if ($value < 0);
  $value = 8192 if ($value > 8192);
  return $value;
}
sub roundoff {
  my ($value) = @_;
  my $rounded_value = int($value + 0.5);
  return $rounded_value;
}

#
# Convert an sRGB value to scRGB, based on IEC 61966-2-2.
#
sub to_scrgb {
  my ($value) = @_;
  my $result;
  if ($value < 10.0) {
    $result = $value * 2.4865;
    return $result;
  }
  $result = ($value + 14.025) / 269.025;
  $result = $result**2.4;
  $result = $result * 8192.0;
  return $result;
}
#
# Convert an scRGB value to sRGB, based on IEC 61966-2-2.
#
sub to_srgb {
  my ($value) = @_;
  my $result;
  if ($value < 0.0) {
    $result = 0.0;
    return $result;
  }
  $result = $value / 8192.0;
  if ($result < 0.00304) {
    $result = 255.0 * (12.92 * $result);
    return $result;
  }
  if ($result <= 1.0) {
    $result = $result ** (1.0 / 2.4);
    $result = $result * 269.025;
    $result = $result - 14.025;
    return $result;
  }
  $result = 255.0;
  return $result;
}

__END__

=head1 NAME

texture_gradient - create a gradient between textures

=head1 SYNOPSIS

texture_gradient [options] output_file

Write a PNG file containing a gradient between textures.

 Options:
  --help        brief help message
  --man         full documentation
  --height      height of the output in pixels, default 2550 
  --width       width of the output in pixels, default 3600
  --verbose     output progress information, can be repeated
  --patch_dir   a directory holding texture patches

=head1 OPTIONS

=over 8

=item B<help>

Print a brief help message and exit.

=item B<man>

Print the manual page and exit.

=item B<height>

The height of the output image, in pixels.  The default is 2550,
which is 8.5 inches at 300 pixels per inch.

=item B<width>

The width of the output image, in pixels.  The default is 3600,
which is 12 inches at 300 pixels per inch.

=item B<verbose>

Output informative and progress messages as the program runs.  This
option can be repeated to gain more information, but messages from
high levels of verbosity are intended for debugging the algorithms,
and can be cryptic.

=item B<patch_dir>

The directory in which to look for patch files, which will be of
the form <name>.png.  The name is a number recording the number of
pixels between the top of the photograph and the center of the patch.

=back

=head1 DESCRIPTION

texture_gradient constructs a gradient between textures.  This is good for 
faking the sky in a photograph.  The anchor points of the gradient are 
constructed from small pieces of the texture that you extract from your photograph.  

Ever had a panorama in which the clouds moved while you were capturing each
image, so in the result the clouds look very unnatural?  Your impulse is to
erase the sky and substitute something for it.  The problem is, sky texture
is very complex, so it is difficult to find a good substiture.  texture_gradient
will construct a good-looking synthetic sky.  Here is how you do it.

Find a column of sky that is mostly uncluttered, and take several small samples
from it.  Note the distance from the top of the photograph, in pixels, of the
center of the sample.  Create a special subdirectory for these samples, and save 
each sample as a PNG (Portable Network Graphics) file, with its file name being
the pixel distance you noted above.  It is OK if the patches are not exactly in
the same column, but be careful not to get two patches of quite different
color with a small vertical distance between them, or you will get a band
between them in the result.  It is best to select clear sky, if you can find it.

Run texture_gradient specifying the height and width of your picture, an
output_file where the texture will be written as a PNG file, and option
patch_dir specifying the directory where you stored the patches.

If you have no sky to match, omit the --patche_dir option and texture_gradient 
will construct a default sky.

Load the output of texture_gradient as an additional layer below your photograph, 
and make the sky transparent.  You should see a clear synthetic sky which matches 
your real sky.  

=head1 BACKGROUND

The color of the sky is remarkably complex.  Each small patch is not a single color, but
a small volume of colors.  That distribution needs to be preserved or the synthesized sky
looks like a cartoon sky.  In addition, the mean color changes with the angle from the
horizon.  All three color components increase in intensity quickly as you begin looking upward from
the horizon, then slowly decrease.  Blue at first levels off at its maximum intensity,
then slowly declines.  Green begins falling more quickly than blue, and falls at a steeper
angle.  Red begins falling more quickly than green, and likewise falls faster at first,
but then begins to level off at high angles.  If the photograph includes both tall and
short objects at a distance, the horizon effect can be observed to follow the contour
of the effective horizon.

=cut