Texture Gradient

Sometimes when I am constructing a panorama, the sky is so messed up due to moving clouds that I just want to replace it. The problem with doing this is that the sky is not a simple color, or even a linear gradient. The color of the sky, even on a clear day, progresses in a complex way from horizon to zenith.

The texture_gradient program will construct a complex gradient from slices that you provide. If you can find good sky material at various heights, it will fill in the transitions.

Even easier, if you don't care about coming close to the original sky, just use the default sky, which matches a late summer day on the Atlantic coast of New Hampshire, as recorded by my Panasonic Lumix DMC-FZ8. Here is that sky:

default sky from texture_gradient

Create a sky to match your picture using --width and --height, then load it as a layer beneath your picture. Erase everything in your picture that should be sky, and the synthesized sky will show instead.

If you want to match your sky and your camera better than the defaults, cut patches of sky and place them in a directory, using file names of the form <position>.png, where "<position>" is the distance in pixels from the top of the picture to the center of the patch. Use --patch_dir to tell texture_gradient where to find the directory.

If you need more sky, texture_gradient will extrapolate. Just add whatever extra height you need to the <position> values, add that much to --height, and increase your picture height by that much before loading the synthesized sky behind it.

texture_gradient also has a --verbose qualifier in case you would like to watch it run. --verbose may be repeated for more logging, but high levels of logging are somewhat cryptic.

The program is written in perl, so should be runnable on GNU/Linux, Microsoft Windows and Apple OS/X, but I have only tested it on GNU/Linux.

Here is the code. You should be able to cut and paste it into a file and execute it from your perl interpreter. If you prefer, you can fetch it from this URL: https://www.systemeyescomputerstore.com/texture_processing/texture_gradient. There is help text at the end of the script.

#!/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