text:
// Implementation of the interface in uint8_image.h.

#include <errno.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <signal.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
#include <assert.h>

#include <glib.h>
#if GLIB_CHECK_VERSION (2, 6, 0)
#  include <glib/gstdio.h>
#endif
#include <gsl/gsl_spline.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_math.h>

#include "asf_jpeg.h"
#include "uint8_image.h"
#include "asf.h"

#ifndef linux
#ifndef darwin
#ifndef win32
static double
round (double arg)
{
    return floor (arg + 0.5);
}
#endif // #ifndef win32
#endif // #ifndef darwin
#endif // #ifndef linux

#include "asf_glib.h"

// Default cache size to use is 16 megabytes.
static const size_t default_cache_size = 16 * 1048576;
// This class wide data element keeps track of the number of temporary
// tile files opened by the current process, in order to give them
// unique names.
static unsigned long current_tile_file_number = 0;
// We need to ensure that multiple threads trying to create their own
// images concurently don't end up with the same temporary file
// names.
//G_LOCK_DEFINE_STATIC (current_tile_file_number);

#ifndef win32
// We don't want to let multiple threads twiddle the signal block mask
// concurrently, or we might end up with the wrong set of signals
// blocked.  This lock is used to guarantee this can't happen (see the
// usage for a better explanation).
G_LOCK_DEFINE_STATIC (signal_block_activity);
#endif

// Return a FILE pointer refering to a new, already unlinked file in a
// location which hopefully has enough free space to serve as a block
// cache.
static FILE *
initialize_tile_cache_file (GString **tile_file_name)
{
  // Create the temporary tile oriented storage file.  This gets
  // filled in in different ways depending on which creation routine
  // we are using.
  g_assert(*tile_file_name == NULL);
  *tile_file_name = g_string_new ("");

  // Here we do a slightly weird thing: if the current directory is
  // writable, we create a temporary file in the current directory.
  // We do this because the temporary file could well be pretty big
  // and /tmp often maps to a small file system.  The idea is that the
  // directory the user is in is more likely to have the extra space
  // required to hold the temporary file.  Of course, if they have
  // been carefully calculating their space requirements, they may be
  // disappointed.  We use a weird name that no sane user would ever
  // use for one of their files, we hope.
  //G_LOCK (current_tile_file_number);
  g_assert (sizeof (long) >= sizeof (pid_t));
  g_string_append_printf (*tile_file_name,
                          ".uint8_image_tile_file_%ld_%lu",
                          (long) getpid (),
                          current_tile_file_number);

  // This hard coded limit on the current number used to uniqueify
  // file names limits us to creating no more than ULONG_MAX instances
  // during a process.
  g_assert (current_tile_file_number < ULONG_MAX);
  current_tile_file_number++;
  //G_UNLOCK (current_tile_file_number);

#ifndef win32
  // We block signals while we create and unlink this file, so we
  // don't end up leaving a huge temporary file somewhere.
  // Theoretically, two parallel instantiations of image could end up
  // in a race condition which would result in all signals ending up
  // blocked after both were done with this section, so we consider
  // this section critical and protect it with a lock.
  G_LOCK (signal_block_activity);
  sigset_t all_signals, old_set;
  int return_code = sigfillset (&all_signals);
  g_assert (return_code == 0);
  return_code = sigprocmask (SIG_SETMASK, &all_signals, &old_set);
#endif

  // FIXME?: It might be faster to use file descriptor based I/O
  // everywhere, or at least for the big transfers.  I'm not sure its
  // worth the trouble though.
  FILE *tile_file = fopen_tmp_file ((*tile_file_name)->str, "w+b");
  if ( tile_file == NULL ) {
    if ( errno != EACCES ) {
      g_warning ("couldn't create file in current directory, and it wasn't"
                 "just a permissions problem");
    }
    else {
      // Couldn't open in current directory, so try using tmpfile,
      // which opens the file in the standardish place for the system.
      // See the comment above about why opening in /tmp or the like
      // is potentially bad.
      tile_file = tmpfile ();
      g_assert (tile_file != NULL);
    }
  }
  else {
#ifndef win32
    return_code = unlink_tmp_file ((*tile_file_name)->str);
    g_assert (return_code == 0);
#endif
  }
  g_assert (tile_file != NULL);

#ifndef win32
  return_code = sigprocmask (SIG_SETMASK, &old_set, NULL);
  G_UNLOCK (signal_block_activity);
#endif

  return tile_file;
}

// This routine does the work common to several of the differenct
// creation routines.  Basically, it does everything but fill in the
// contents of the disk tile store.
static UInt8Image *
initialize_uint8_image_structure (ssize_t size_x, ssize_t size_y)
{
  // Allocate instance memory.
  UInt8Image *self = g_new0 (UInt8Image, 1);

  // Validate and remember image size.
  g_assert (size_x > 0 && size_y > 0);
  self->size_x = size_x;
  self->size_y = size_y;

  // Greater of size_x and size_y.
  size_t largest_dimension = (size_x > size_y ? size_x : size_y);

  // If we can fit the entire image in a single square tile, then we
  // want just a single big tile and we won't need to bother with the
  // cache file since it won't ever be used, so we do things slightly
  // differently.  FIXME: it would be slightly better to also detect
  // and specially handle the case where we have long narrow images
  // that can fit in a single stip of tiles in the cache.
  if ( largest_dimension * largest_dimension * sizeof (uint8_t)
       <= default_cache_size ) {
    self->cache_space = (largest_dimension * largest_dimension
       * sizeof (uint8_t));
    self->cache_area = self->cache_space / sizeof (uint8_t);
    self->tile_size = largest_dimension;
    self->cache_size_in_tiles = 1;
    self->tile_count_x = 1;
    self->tile_count_y = 1;
    self->tile_count = 1;
    self->tile_area = self->tile_size * self->tile_size;
    self->cache = g_new (uint8_t, self->cache_area);
    self->tile_addresses = g_new0 (uint8_t *, self->tile_count);
    g_assert (NULL == 0x0);     // Ensure g_new0 effectively sets to NULL.
    // The tile queue shouldn't ever be needed in this case.
    self->tile_queue = NULL;
    // The tile file shouldn't ever be needed, so we set it to NULL to
    // indicate this to a few other methods that use it directly, and
    // to hopefully ensure that it triggers an exception if it is
    // used.
    self->tile_file = NULL;

    return self;
  }

  // The default cache size compiled into the class.
  self->cache_space = default_cache_size;

  // Memory cache space, in pixels.
  g_assert (self->cache_space % sizeof (uint8_t) == 0);
  self->cache_area = self->cache_space / sizeof (uint8_t);

  // How small do our tiles have to be on a side to fit two full rows
  // of them in the memory cache?  This is slightly tricky.  In order
  // to provide the services promised in the interface, we need to
  // solve
  //
  //      2 * pow (t, 2) * ceil ((double)largest_dimension / t)
  //           <= self->cache_area
  //
  // for tile size t.  I don't know the closed form solution if there
  // is one, so toss out the ceil() and solve the easier
  //
  //      2 * pow (t, 2) * ((double)largest_dimension / t) <= self->cache_area
  //
  // and then decrement t iteratively until things work.
  self->tile_size = self->cache_area / (2 * largest_dimension);
  while ( (2 * self->tile_size * self->tile_size
           * ceil ((double) largest_dimension / self->tile_size))
          > self->cache_area ) {
    self->tile_size--;
  }

  // Area of tiles, in pixels.
  self->tile_area = (size_t) (self->tile_size * self->tile_size);

  // Number of tiles which will fit in image cache.
  self->cache_size_in_tiles = self->cache_area / self->tile_area;
  // Can we fit at least as much as we intended in the cache?
  g_assert (self->cache_size_in_tiles
            >= 2 * (size_t) ceil ((double) largest_dimension
                                  / self->tile_size));

  // Number of tiles image has been split into in x and y directions.
  self->tile_count_x = (size_t) ceil ((double) self->size_x / self->tile_size);
  self->tile_count_y = (size_t) ceil ((double) self->size_y / self->tile_size);

  // Total number of image tiles image has been split into.
  self->tile_count = self->tile_count_x * self->tile_count_y;

  // We want to be able to pack a tile number into a pointer later, so
  // we need it to fit into an integer.
  g_assert (self->tile_count < INT_MAX);

  // Did all that math work?
  g_assert (self->tile_size * self->cache_size_in_tiles / 2
            >= largest_dimension);
  g_assert (self->cache_size_in_tiles * self->tile_area <= self->cache_area);

  // Allocate memory for the in-memory cache.
  self->cache = g_new (uint8_t, self->cache_area);
  // Do we want to do mlock() here maybe?

  // The addresses in the cache of the starts of each of the tiles.
  // This array contains flattened tile addresses in the same way that
  // image memory normally uses flattened pixel addresses, e.g. the
  // address of tile x = 2, y = 4 is stored at self->tile_addresses[4
  // * self->tile_count_x + 2].  If a tile isn't in the cache, the
  // address is NULL (meaning it will have to be loaded).
  self->tile_addresses = g_new0 (uint8_t *, self->tile_count);
  g_assert (NULL == 0x0);       // Ensure g_new0 effectively sets to NULL.

  // Create a queue in order to keep track of which tile was loaded
  // longest ago.
  self->tile_queue = g_queue_new ();

  // Get a new empty tile cache file pointer.
  self->tile_file_name = NULL;
  self->tile_file = initialize_tile_cache_file ( &(self->tile_file_name) );

  return self;
}

UInt8Image *
uint8_image_thaw (FILE *file_pointer)
{
  FILE *fp = file_pointer;  // Convenience alias.

  g_assert (file_pointer != NULL);

  UInt8Image *self = g_new (UInt8Image, 1);

  size_t read_count = fread (&(self->size_x), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->size_y), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->cache_space), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->cache_area), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->tile_size), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->cache_size_in_tiles), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->tile_count_x), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->tile_count_y), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->tile_count), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  read_count = fread (&(self->tile_area), sizeof (size_t), 1, fp);
  g_assert (read_count == 1);

  // The cache isn't serialized -- its a bit of a pain and probably
  // almost never worth it.
  self->cache = g_new (uint8_t, self->cache_area);

  self->tile_addresses = g_new0 (uint8_t *, self->tile_count);

  // We don't actually keep the tile queue in the serialized instance,
  // but if the serialized pointer to it is NULL, we know we aren't
  // using a tile cache file (i.e. the whole image fits in the memory
  // cache).
  read_count = fread (&(self->tile_queue), sizeof (GQueue *), 1, fp);
  g_assert (read_count == 1);

  // If there was no cache file...
  if ( self->tile_queue == NULL ) {
    // The tile_file structure field should also be NULL.
    self->tile_file = NULL;
    // we restore the file directly into the first and only tile (see
    // the end of the uint8_image_new method).
    self->tile_addresses[0] = self->cache;
    read_count = fread (self->tile_addresses[0], sizeof (uint8_t),
      self->tile_area, fp);
    g_assert (read_count == self->tile_area);
  }
  // otherwise, an empty tile queue needs to be initialized, and the
  // remainder of the serialized version is the tile block cache.
  else {
    self->tile_queue = g_queue_new ();
    self->tile_file_name = NULL;
    self->tile_file = initialize_tile_cache_file ( &(self->tile_file_name) );
    uint8_t *buffer = g_new (uint8_t, self->tile_area);
    size_t ii;
    for ( ii = 0 ; ii < self->tile_count ; ii++ ) {
      read_count = fread (buffer, sizeof (uint8_t), self->tile_area, fp);
      g_assert (read_count == self->tile_area);
      size_t write_count = fwrite (buffer, sizeof (uint8_t), self->tile_area,
           self->tile_file);
      if ( write_count < self->tile_area ) {
  if ( feof (self->tile_file) ) {
    fprintf (stderr,
       "Premature end of file while trying to thaw UInt8Image "
       "instance\n");
  }
  else {
    g_assert (ferror (self->tile_file));
    fprintf (stderr,
       "Error writing tile cache file for UInt8Image instance "
       "during thaw: %s\n", strerror (errno));
  }
  exit (EXIT_FAILURE);
      }
      g_assert (write_count == self->tile_area);
    }
    g_free (buffer);
  }

  return self;
}

UInt8Image *
uint8_image_new (ssize_t size_x, ssize_t size_y)
{
  g_assert (size_x > 0 && size_y > 0);

  UInt8Image *self = initialize_uint8_image_structure (size_x, size_y);

  // If we need a tile file for an image of this size, prepare it.
  if ( self->tile_file != NULL ) {

    // The total width or height of all the tiles is probably greater
    // than the width or height of the image itself.
    size_t total_width = self->tile_count_x * self->tile_size;
    size_t total_height = self->tile_count_y * self->tile_size;

    // Fill the file full of zeros.  FIXME: there is almost certainly
    // a faster way to ensure that we have the disk space we need.
    uint8_t *zero_line = g_new0 (uint8_t, total_width);
    g_assert (0 == 0x0);      // Ensure the g_new0 did what we think.

    // We don't have to write in tile order because its all zeros anyway.
    size_t ii;
    for ( ii = 0 ; ii < total_height ; ii++ ) {
      size_t write_count = fwrite (zero_line, sizeof (uint8_t), total_width,
                                   self->tile_file);
      // If we wrote less than expected,
      if ( write_count < total_width ) {
        // it must have been a write error (probably no space left),
        g_assert (ferror (self->tile_file));
        // so print an error message,
        fprintf (stderr,
                 "Error writing tile cache file for UInt8Image instance: "
                 "%s\n", strerror (errno));
        // and exit.
        exit (EXIT_FAILURE);
      }
    }

    // Done with the line of zeros.
    g_free (zero_line);
  }

  // Everything fits in the cache (at the moment this means everything
  // fits in the first tile, which is a bit of a FIXME), so just put
  // it there.
  else {
    self->tile_addresses[0] = self->cache;
    size_t ii, jj;
    for ( ii = 0 ; ii < self->tile_size ; ii++ ) {
      for ( jj = 0 ; jj < self->tile_size ; jj++ ) {
        self->tile_addresses[0][ii * self->tile_size + jj] = 0;
      }
    }
  }

  return self;
}

UInt8Image *
uint_image_new_with_value (ssize_t size_x, ssize_t size_y, uint8_t value)
{
  // Carefully clone-and-modified over from float_image.c, but not
  // tested yet.
  g_assert_not_reached ();

  g_assert (size_x > 0 && size_y > 0);

  UInt8Image *self = initialize_uint8_image_structure (size_x, size_y);

  // If we need a tile file for an image of this size, prepare it.
  if ( self->tile_file != NULL ) {

    // The total width or height of all the tiles is probably greater
    // than the width or height of the image itself.
    size_t total_width = self->tile_count_x * self->tile_size;
    size_t total_height = self->tile_count_y * self->tile_size;

    // Fill the file full of the given value.
    uint8_t *value_line = g_new (uint8_t, total_width);
    size_t ii;
    for ( ii = 0 ; ii < total_width ; ii++ ) {
      value_line[ii] = value;
    }
    // We don't have to write in tile order because the values are all
    // the same anyway.
    for ( ii = 0 ; ii < total_height ; ii++ ) {
      size_t write_count = fwrite (value_line, sizeof (uint8_t), total_width,
                                   self->tile_file);
      // If we wrote less than expected,
      if ( write_count < total_width ) {
        // it must have been a write error (probably no space left),
        g_assert (ferror (self->tile_file));
        // so print an error message,
        fprintf (stderr,
                 "Error writing tile cache file for UInt8Image instance: "
                 "%s\n", strerror (errno));
        // and exit.
        exit (EXIT_FAILURE);
      }
    }

    // Done with the line of values.
    g_free (value_line);
  }

  // Everything fits in the cache (at the moment this means everything
  // fits in the first tile, which is a bit of a FIXME), so just put
  // it there.
  else {
    self->tile_addresses[0] = self->cache;
    size_t ii, jj;
    for ( ii = 0 ; ii < self->tile_size ; ii++ ) {
      for ( jj = 0 ; jj < self->tile_size ; jj++ ) {
        self->tile_addresses[0][ii * self->tile_size + jj] = value;
      }
    }
  }

  return self;
}

UInt8Image *
uint8_image_new_from_memory (ssize_t size_x, ssize_t size_y, uint8_t *buffer)
{
  g_assert (size_x > 0 && size_y > 0);

  g_assert_not_reached ();      // Stubbed out for now.
  // Compiler reassurance.
  size_x = size_x;
  size_y = size_y;
  buffer = buffer;
  return NULL;
}

// Bilinear interpolation for a point delta_x, delta_y from the lower
// left corner between values ul (upper left), ur (upper right), etc.
// The corner are considered to be corners of a unit square.
static float
bilinear_interpolate (double delta_x, double delta_y, float ul, float ur,
                      float ll, float lr)
{
  float lv = ll + (lr - ll) * delta_x; // Lower value.
  float uv = ul + (ur - ul) * delta_x; // Upper value.

  return lv + (uv - lv) * delta_y;
}

UInt8Image *
uint8_image_copy (UInt8Image *model)
{
  // This method should be totally fine, it is extremely trivial clone
  // and modify from the corresponding method in float_image.c, but
  // since its untested I disable it for the moment.
  g_assert_not_reached ();

  // FIXME: this could obviously be optimized a lot by copying the
  // existed tile file, etc.
  UInt8Image *self = uint8_image_new (model->size_x, model->size_y);

  size_t ii, jj;
  for ( ii = 0 ; ii < self->size_y ; ii++ ) {
    for ( jj = 0 ; jj < self->size_x ; jj++ ) {
      uint8_image_set_pixel (self, jj, ii,
           uint8_image_get_pixel (model, jj, ii));
    }
  }

  return self;
}

UInt8Image *
uint8_image_new_from_model_scaled (UInt8Image *model, ssize_t scale_factor)
{
  g_assert (model->size_x > 0 && model->size_y > 0);

  g_assert (scale_factor > 0);
  g_assert (scale_factor % 2 == 1);

  UInt8Image *self
    = uint8_image_new (round ((double) model->size_x / scale_factor),
           round ((double) model->size_y / scale_factor));

  // This method hasn't yet made the jump from FloatImage to here.
  // The implementation should ultimately follow the one in FloatImage
  // to achieve the interface description given in uint8_image.h.
  g_assert_not_reached ();

  return self;
}

UInt8Image *
uint8_image_new_subimage (UInt8Image *model, ssize_t x, ssize_t y,
        ssize_t size_x, ssize_t size_y)
{
  // Upper left corner must be in model.
  g_assert (x >= 0 && y >= 0);

  // Size of image to be created must be strictly positive.
  g_assert (size_x >= 1 && size_y >= 1);

  // Given model must be big enough to allow a subimage of the
  // requested size to fit.
  g_assert (model->size_x <= SSIZE_MAX && model->size_y <= SSIZE_MAX);
  g_assert (x + size_x <= (ssize_t) model->size_x);
  g_assert (y + size_y <= (ssize_t) model->size_y);

  UInt8Image *self = uint8_image_new (size_x, size_y);

  // Copy the image pixels from the model.
  ssize_t ii, jj;
  for ( ii = 0 ; ii < (ssize_t) self->size_x ; ii++ ) {
    for ( jj = 0 ; jj < (ssize_t) self->size_y ; jj++ ) {
      uint8_t pv = uint8_image_get_pixel (model, x + ii, y + jj);
      uint8_image_set_pixel (self, ii, jj, pv);
    }
  }

  return self;
}

// Return true iff file is size or larger.
static gboolean
is_large_enough (const char *file, off_t size)
{
  struct stat stat_buffer;
#if GLIB_CHECK_VERSION(2, 6, 0)
  int return_code = g_stat (file, &stat_buffer);
  if ( return_code != 0 ) {
    g_error ("Couldn't g_stat file %s: %s", file, strerror (errno));
  }
#else
  int return_code = stat (file, &stat_buffer);
  if ( return_code != 0 ) {
    g_error ("Couldn't stat file %s: %s", file, strerror (errno));
  }
#endif

  return stat_buffer.st_size >= size;
}

UInt8Image *
uint8_image_new_from_file (ssize_t size_x, ssize_t size_y, const char *file,
                           off_t offset)
{
  g_assert (size_x > 0 && size_y > 0);

  // Check in advance if the source file looks big enough (we will
  // still need to check return codes as we read() data, of course).
  g_assert (is_large_enough (file, offset + ((off_t) size_x * size_y
               * sizeof (uint8_t))));

  // Open the file to read data from.
  FILE *fp = fopen (file, "rb");
  // FIXME: we need some error handling and propagation here.
  g_assert (fp != NULL);

  UInt8Image *self = uint8_image_new_from_file_pointer (size_x, size_y, fp,
                                                        offset);

  // Close file we read image from.
  int return_code = fclose (fp);
  g_assert (return_code == 0);

  return self;
}

// Return true iff file referred to by file_pointer is larger than size.
static gboolean
file_pointed_to_larger_than (FILE *file_pointer, off_t size)
{
  struct stat stat_buffer;
  int return_code = fstat (fileno (file_pointer), &stat_buffer);
  g_assert (return_code == 0);
  return stat_buffer.st_size >= size;
}

UInt8Image *
uint8_image_new_from_file_pointer (ssize_t size_x, ssize_t size_y,
                                   FILE *file_pointer, off_t offset)
{
  g_assert (size_x > 0 && size_y > 0);

  // Check in advance if the source file looks big enough (we will
  // still need to check return codes as we read() data, of course).
  g_assert (file_pointed_to_larger_than (file_pointer,
           offset + ((off_t) size_x * size_y
               * sizeof (uint8_t))));

  UInt8Image *self = initialize_uint8_image_structure (size_x, size_y);

  FILE *fp = file_pointer;      // Convenience alias.

  // Seek to the indicated offset in the file.
  int return_code = FSEEK64 (fp, offset, SEEK_CUR);
  g_assert (return_code == 0);

  // If we need a tile file for an image of this size, we will load
  // the data straight into it.
  if ( self->tile_file != NULL ) {

    // We will read the input image data in horizontal stips one tile
    // high.  Note that we probably won't be able to entirely fill the
    // last tiles in each dimension with real data, since the image
    // sizes rarely divide evenly by the numbers of tiles.  So we fill
    // it with zeros instead.  The data off the edges of the image
    // should never be accessed directly anyway.

    // Some data for doing zero fill.  If the tiles are bigger than
    // the image itself, we need to make the available zero fill the
    // size of the tile instead of the size of the file.
    g_assert (self->tile_size <= SSIZE_MAX);
    uint8_t *zero_line = g_new0 (uint8_t, (size_x > (ssize_t) self->tile_size ?
             (size_t) size_x : self->tile_size));
    g_assert (0 == 0x0);

    // Buffer capable of holding a full strip.
    uint8_t *buffer = g_new (uint8_t, self->tile_size * self->size_x);

    // Reorganize data into tiles in tile oriented disk file.
    size_t ii = 0;
    for ( ii = 0 ; ii < self->tile_count_y ; ii++ ) {
      // The "effective_height" of the strip is the portion of the
      // strip for which data actually exists.  If the effective
      // height is less than self->tile>size, we will have to add some
      // junk to fill up the extra part of the tile (which should
      // never be accessed).
      size_t effective_height;
      if ( ii < self->tile_count_y - 1
           || self->size_y % self->tile_size == 0 ) {
        effective_height = self->tile_size;
      }
      else {
        effective_height = self->size_y % self->tile_size;
      }
      // Total area of the current strip.
      size_t strip_area = effective_height * self->size_x;

      // Read one strip of tiles worth of data from the file.
      size_t read_count = fread (buffer, sizeof (uint8_t), strip_area, fp);
      g_assert (read_count == strip_area);

      // Write data from the strip into the tile store.
      size_t jj;
      for ( jj = 0 ; jj < self->tile_count_x ; jj++ ) {
        // This is roughly analogous to effective_height.
        size_t effective_width;
        if ( jj < self->tile_count_x - 1
             || self->size_x % self->tile_size == 0) {
          effective_width = self->tile_size;
        }
        else {
          effective_width = self->size_x % self->tile_size;
        }
        size_t write_count;     // For return of fwrite() calls.
        size_t kk;
        for ( kk = 0 ; kk < effective_height ; kk++ ) {
          write_count
            = fwrite (buffer + kk * self->size_x + jj * self->tile_size,
                      sizeof (uint8_t), effective_width, self->tile_file);
          // If we wrote less than expected,
          if ( write_count < effective_width ) {
            // it must have been a write error (probably no space left),
            g_assert (ferror (self->tile_file));
            // so print an error message,
            fprintf (stderr,
                     "Error writing tile cache file for UInt8Image "
                     "instance: %s\n", strerror (errno));
            // and exit.
            exit (EXIT_FAILURE);
          }
          if ( effective_width < self->tile_size ) {
            // Amount we have left to write to fill out the last tile.
            size_t edge_width = self->tile_size - effective_width;
            write_count = fwrite (zero_line, sizeof (uint8_t), edge_width,
                                  self->tile_file);
            // If we wrote less than expected,
            if ( write_count < edge_width ) {
              // it must have been a write error (probably no space left),
              g_assert (ferror (self->tile_file));
              // so print an error message,
              fprintf (stderr,
                       "Error writing tile cache file for UInt8Image "
                       "instance: %s\n", strerror (errno));
              // and exit.
              exit (EXIT_FAILURE);
            }
          }
        }
        // Finish writing the bottom of the tile for which there is no
        // image data (should only happen if we are on the last strip of
        // tiles).
        for ( ; kk < self->tile_size ; kk++ ) {
          g_assert (ii == self->tile_count_y - 1);
          write_count = fwrite (zero_line, sizeof (uint8_t), self->tile_size,
                                self->tile_file);
          // If we wrote less than expected,
          if ( write_count < self->tile_size ) {
            // it must have been a write error (probably no space left),
            g_assert (ferror (self->tile_file));
            // so print an error message,
            fprintf (stderr,
                     "Error writing tile cache file for UInt8Image "
                     "instance: %s\n", strerror (errno));
            // and exit.
            exit (EXIT_FAILURE);
          }
        }
      }
    }

    // Did we write the correct total amount of data?
    g_assert (FTELL64 (self->tile_file)
        == ((off_t) self->tile_area * self->tile_count
      * sizeof (uint8_t)));

    // Free temporary buffers.
    g_free (buffer);
    g_free (zero_line);
  }

  // Everything fits in the cache (at the moment this means everything
  // fits in the first tile, which is a bit of a FIXME), so just put
  // it there.
  else {
    self->tile_addresses[0] = self->cache;

    size_t ii;
    for ( ii = 0 ; ii < self->size_y ; ii++ ) {
      // Address where the current row of pixels should end up.
      uint8_t *row_address = self->tile_addresses[0] + ii * self->tile_size;

      // Read the data.
      size_t read_count = fread (row_address, sizeof (uint8_t), self->size_x,
         fp);
      g_assert (read_count == self->size_x);
    }
  }

  return self;
}

UInt8Image *
uint8_image_new_from_file_scaled (ssize_t size_x, ssize_t size_y,
                                  ssize_t original_size_x,
                                  ssize_t original_size_y,
          const char *file, off_t offset)
{
  // This method has been carefully translated from the corresponding
  // method in float_image.c, but it hasn't been tested yet.
  g_assert_not_reached ();

  g_assert (size_x > 0 && size_y > 0);
  g_assert (original_size_x > 0 && original_size_y > 0);

  // Image can only be scaled down with this routine, not up.
  g_assert (size_x < original_size_x);
  g_assert (size_y < original_size_y);

  // Check in advance if the source file looks big enough (we will
  // still need to check return codes as we read() data, of course).
  g_assert (is_large_enough (file,
           offset + ((off_t) original_size_x
               * original_size_y
               * sizeof (uint8_t))));

  // Find the stride that we need to use in each dimension to evenly
  // cover the original image space.
  double stride_x = (size_x == 1 ? 0.0
         : (double) (original_size_x - 1) / (size_x - 1));
  double stride_y = (size_y == 1 ? 0.0
         : (double) (original_size_y - 1) / (size_y - 1));

  // Open the file to read data from.
  FILE *fp = fopen (file, "rb");
  // FIXME: we need some error handling and propagation here.
  g_assert (fp != NULL);

  // We will do a row at a time to save some possibly expensive
  // seeking.  So here we have an entire row worth of upper lefts,
  // upper rights, etc.
  uint8_t *uls = g_new (uint8_t, size_x);
  uint8_t *urs = g_new (uint8_t, size_x);
  uint8_t *lls = g_new (uint8_t, size_x);
  uint8_t *lrs = g_new (uint8_t, size_x);

  // Results of rounded bilinear interpolation for the current row.
  uint8_t *interpolated_values = g_new (uint8_t, size_x);

  // We will write the reduced resolution version of the image into a
  // temporary file so we can leverage the new_from_file method and
  // avoid trying to stick the whole reduced resolution image in
  // memory.
  FILE *reduced_image = tmpfile ();

  ssize_t ii, jj;
  for ( ii = 0 ; ii < size_y ; ii++ ) {
    size_t read_count;          // For fread calls.
    int return_code;            // For FSEEK64 calls.
    // Input image y index of row above row of interest.
    ssize_t in_ray = floor (ii * stride_y);
    // Due to the vagaries of floating point arithmetic, we might run
    // past the index of our last pixel by a little bit, so we correct.
    if ( in_ray >= original_size_y - 1 ) {
      // We better not be much over the last index though.
      g_assert (in_ray < original_size_y);
      // The index should be an integer, so floor should fix us up.
      in_ray = floor (in_ray);
      g_assert (in_ray == original_size_y - 1);
    }
    g_assert (in_ray < original_size_y);
    // Input image y index of row below row of interest.  If we would
    // be off the image, we just take the last row a second time, and
    // let the interpolation work things out.
    ssize_t in_rby;
    if ( in_ray == original_size_y - 1 ) {
      in_rby = in_ray;
    }
    else {
      in_rby = in_ray + 1;
    }
    // Fetch the row above.
    for ( jj = 0 ; jj < size_x ; jj++ ) {
      // Input image indicies of current upper left corner pixel.
      ssize_t in_ul_x = floor (jj * stride_x);
      // Watch for floating point inexactness (see comment above).
      if ( G_UNLIKELY (in_ul_x >= original_size_x - 1) ) {
        g_assert (in_ul_x < original_size_x);
        in_ul_x = floor (in_ul_x);
        g_assert (in_ul_x == original_size_x - 1);
      }
      g_assert (in_ul_x < original_size_x);
      size_t in_ul_y = in_ray;
      off_t sample_offset
  = offset + sizeof (uint8_t) * ((off_t) in_ul_y * original_size_x
               + in_ul_x);
      return_code = FSEEK64 (fp, sample_offset, SEEK_SET);
      g_assert (return_code == 0);
      read_count = fread (&(uls[jj]), sizeof (uint8_t), 1, fp);
      g_assert (read_count == 1);
      // If the upper left pixel was the last pixel in the input image,
      if ( in_ul_x == original_size_x - 1 ) {
        // just treat it as the upper right as well,
        urs[jj] = uls[jj];
      }
      // otherwise read the next pixel as the upper right pixel.
      else {
        read_count = fread (&(urs[jj]), sizeof (uint8_t), 1, fp);
        g_assert (read_count == 1);
      }
    }
    // Fetch the row below.
    for ( jj = 0 ; jj < size_x ; jj++ ) {
      // Input image indicies of the lower left corner pixel.
      ssize_t in_ll_x = floor (jj * stride_x);
      // Watch for floating point inexactness (see comment above).
      if ( G_UNLIKELY (in_ll_x >= original_size_y - 1) ) {
        g_assert (in_ll_x < original_size_x);
        in_ll_x = floor (in_ll_x);
        g_assert (in_ll_x == original_size_x - 1);
      }
      g_assert (in_ll_x < original_size_x);
      size_t in_ll_y = in_rby;
      off_t sample_offset
  = offset + sizeof (uint8_t) * ((off_t) in_ll_y * original_size_x
               + in_ll_x);
      return_code = FSEEK64 (fp, sample_offset, SEEK_SET);
      g_assert (return_code == 0);
      read_count = fread (&(lls[jj]), sizeof (uint8_t), 1, fp);
      g_assert (read_count == 1);
      // If the lower left pixel was the last pixel in the input image,
      if ( in_ll_x == original_size_x - 1 ) {
        // just treat it as the lower right as well,
        lrs[jj] = lls[jj];
      }
      // otherwise read the next pixel as the lower right pixel.
      else {
        read_count = fread (&(lrs[jj]), sizeof (uint8_t), 1, fp);
        g_assert (read_count == 1);
      }
    }

    // Perform the interpolation.
    for ( jj = 0 ; jj < size_x ; jj++ ) {
      double delta_x = stride_x * jj - floor (stride_x * jj);
      double delta_y = -(stride_y * ii - floor (stride_y * ii));
      float interpolated_value = bilinear_interpolate (delta_x, delta_y,
                   uls[jj], urs[jj],
                   lls[jj], lrs[jj]);
      g_assert (interpolated_value >= 0);
      g_assert (interpolated_value <= UINT8_MAX);
      interpolated_values[jj] = round (interpolated_value);
    }
    size_t write_count = fwrite (interpolated_values, sizeof (uint8_t), size_x,
                                 reduced_image);
    g_assert (write_count == (size_t) size_x);
  }

  // We are done with the temporary buffers.
  g_free (interpolated_values);
  g_free (lrs);
  g_free (lls);
  g_free (urs);
  g_free (uls);

  // Reposition to the beginning of the temporary file to fit with
  // operation of new_from_file_pointer method.
  int return_code = FSEEK64 (reduced_image, (off_t) 0, SEEK_SET);
  g_assert (return_code == 0);

  // Slurp the scaled file back in as an instance.
  UInt8Image *self
    = uint8_image_new_from_file_pointer (size_x, size_y, reduced_image,
           (off_t) 0);

  // Now that we have an instantiated version of the image we are done
  // with this temporary file.
  return_code = fclose (reduced_image);

  return self;
}

// Returns a new UInt8Image, for the image corresponding to the given metadata.
UInt8Image *
uint8_image_new_from_metadata(meta_parameters *meta, const char *file)
{
  return uint8_image_band_new_from_metadata(meta, 0, file);
}

// Returns a new UInt8Image, for the image band corresponding to the
// given metadata.
UInt8Image *
uint8_image_band_new_from_metadata(meta_parameters *meta,
           int band, const char *file)
{
    int nl = meta->general->line_count;
    int ns = meta->general->sample_count;

    FILE * fp = FOPEN(file, "rb");
    UInt8Image * bi = uint8_image_new(ns, nl);

    int i,j;
    unsigned char *buf = MALLOC(sizeof(unsigned char)*ns);
    for (i = 0; i < nl; ++i) {
        get_byte_line(fp, meta, i+band*nl, buf);
        for (j = 0; j < ns; ++j)
            uint8_image_set_pixel(bi, j, i, buf[j]);
    }

    free(buf);
    fclose(fp);

    return bi;
}

// Copy the contents of tile with flattened offset tile_offset from
// the memory cache to the disk file.  Its probably easiest to
// understand this function by looking at how its used.
static void
cached_tile_to_disk (UInt8Image *self, size_t tile_offset)
{
  // If we aren't using a tile file, this operation doesn't make
  // sense.
  g_assert (self->tile_file != NULL);

  // We must have a legitimate tile_offset.
  g_assert (tile_offset < self->tile_count);

  // The tile we are trying to copy from cache to disk must be loaded
  // in the cache for this operation to make sense.
  g_assert (self->tile_addresses[tile_offset] != NULL);

  int return_code
    = FSEEK64 (self->tile_file,
        (off_t) tile_offset * self->tile_area * sizeof (uint8_t),
        SEEK_SET);
  g_assert (return_code == 0);
  size_t write_count = fwrite (self->tile_addresses[tile_offset],
             sizeof (uint8_t), self->tile_area,
             self->tile_file);
  g_assert (write_count == self->tile_area);
}

// Return true iff tile (x, y) is already loaded into the memory cache.
static gboolean
tile_is_loaded (UInt8Image *self, ssize_t x, ssize_t y)
{
  g_assert (    x >= 0 && (size_t) x < self->tile_count_x
             && y >= 0 && (size_t) y < self->tile_count_y );

  size_t tile_offset = self->tile_count_x * y + x;

  return self->tile_addresses[tile_offset] != NULL;
}

// Load (currently unloaded) tile (x, y) from disk cache into memory
// cache, possibly displacing the oldest tile already loaded, updating
// the load order queue, and returning the address of the tile loaded.
static uint8_t *
load_tile (UInt8Image *self, ssize_t x, ssize_t y)
{
  // Make sure we haven't screwed up somehow and not created a tile
  // file when in fact we should have.
  g_assert (self->tile_file != NULL);

  g_assert (!tile_is_loaded (self, x, y));

  // Address into which tile gets loaded (to be returned).
  uint8_t *tile_address;

  // Offset of tile in flattened array.
  size_t tile_offset = self->tile_count_x * y + x;

  // We have to check and see if we have to displace an already loaded
  // tile or not.
  if ( self->tile_queue->length == self->cache_size_in_tiles ) {
    // Displace tile loaded longest ago.
    size_t oldest_tile
      = GPOINTER_TO_INT (g_queue_pop_tail (self->tile_queue));
    cached_tile_to_disk (self, oldest_tile);
    tile_address = self->tile_addresses[oldest_tile];
    self->tile_addresses[oldest_tile] = NULL;
  }
  else {
    // Load tile into first free slot.
    tile_address = self->cache + self->tile_queue->length * self->tile_area;
  }

  // Put the new tile address into the index, and put the index into
  // the load order queue.
  self->tile_addresses[tile_offset] = tile_address;
  // Stash in queue by converting to a pointer (so it must fit in an int).
  g_assert (tile_offset < INT_MAX);
  g_queue_push_head (self->tile_queue,
                     GINT_TO_POINTER ((int) tile_offset));

  // Load the tile data.
  int return_code
    = FSEEK64 (self->tile_file,
              (off_t) tile_offset * self->tile_area * sizeof (uint8_t),
              SEEK_SET);
  g_assert (return_code == 0);
  clearerr (self->tile_file);
  size_t read_count = fread (tile_address, sizeof (uint8_t), self->tile_area,
                             self->tile_file);
  if ( read_count < self->tile_area ) {
    if ( ferror (self->tile_file) ) {
      perror ("error reading tile cache file");
      g_assert_not_reached ();
    }
    if ( feof (self->tile_file) ) {
      fprintf (stderr,
               "nothing left to read in tile cache file at offset %lld\n",
               FTELL64 (self->tile_file));
      g_assert_not_reached ();
    }
  }
  g_assert (read_count == self->tile_area);

  return tile_address;
}

uint8_t
uint8_image_get_pixel (UInt8Image *self, ssize_t x, ssize_t y)
{
  // Are we at a valid image pixel?
  g_assert (x >= 0 && (size_t) x < self->size_x);
  g_assert (y >= 0 && (size_t) y < self->size_y);

  // Get the pixel coordinates, including tile and pixel-in-tile.
  g_assert (sizeof (long int) >= sizeof (size_t));
  ldiv_t pc_x = ldiv (x, self->tile_size), pc_y = ldiv (y, self->tile_size);

  // Offset of tile x, y, where tiles are viewed as pixels normally are.
  size_t tile_offset = self->tile_count_x * pc_y.quot + pc_x.quot;

  // Address of data for tile containing pixel of interest (may still
  // have to be loaded from disk cache).
  uint8_t *tile_address = self->tile_addresses[tile_offset];

  // Load the tile containing the pixel of interest if necessary.
  if ( G_UNLIKELY (tile_address == NULL) ) {
    tile_address = load_tile (self, pc_x.quot, pc_y.quot);
  }

  // Return pixel of interest.
  return tile_address[self->tile_size * pc_y.rem + pc_x.rem];
}

void
uint8_image_set_pixel (UInt8Image *self, ssize_t x, ssize_t y, uint8_t value)
{
  // Are we at a valid image pixel?
  g_assert (self != NULL);
  g_assert (x >= 0 && (size_t) x <= self->size_x);
  g_assert (y >= 0 && (size_t) y <= self->size_y);

  // Get the pixel coordinates, including tile and pixel-in-tile.
  g_assert (sizeof (long int) >= sizeof (size_t));
  ldiv_t pc_x = ldiv (x, self->tile_size), pc_y = ldiv (y, self->tile_size);

  // Offset of tile x, y, where tiles are viewed as pixels normally are.
  size_t tile_offset = self->tile_count_x * pc_y.quot + pc_x.quot;

  // Address of data for tile containing pixel of interest (may still
  // have to be loaded from disk cache).
  uint8_t *tile_address = self->tile_addresses[tile_offset];

  // Load the tile containing the pixel of interest if necessary.
  if ( G_UNLIKELY (tile_address == NULL) ) {
    tile_address = load_tile (self, pc_x.quot, pc_y.quot);
  }

  // Set pixel of interest.
  tile_address[self->tile_size * pc_y.rem + pc_x.rem] = value;
}

void
uint8_image_get_region (UInt8Image *self, ssize_t x, ssize_t y, ssize_t size_x,
                        ssize_t size_y, uint8_t *buffer)
{
  g_assert (size_x >= 0);
  g_assert (x >= 0);
  g_assert ((size_t) x + (size_t) size_x - 1 < self->size_x);
  g_assert (size_y >= 0);
  g_assert (y >= 0);
  g_assert ((size_t) y + (size_t) size_y - 1 < self->size_y);

  ssize_t ii, jj;               // Index variables.
  for ( ii = 0 ; ii < size_y ; ii++ ) {
    for ( jj = 0 ; jj < size_x ; jj++ ) {
      // We are essentially returning a subimage from the big image.
      // These are the indicies in the big image (self) of the current
      // pixel.
      size_t ix = x + jj, iy = y + ii;
      buffer[ii * size_x + jj] = uint8_image_get_pixel (self, ix, iy);
    }
  }
}

void
uint8_image_set_region (UInt8Image *self, size_t x, size_t y, size_t size_x,
                        size_t size_y, uint8_t *buffer)
{
  g_assert_not_reached ();      // Stubbed out for now.
  self = self; x = x; y = y; size_x = size_x, size_y = size_y; buffer = buffer;
}

void
uint8_image_get_row (UInt8Image *self, size_t row, uint8_t *buffer)
{
  uint8_image_get_region (self, 0, row, self->size_x, 1, buffer);
}

uint8_t
uint8_image_get_pixel_with_reflection (UInt8Image *self, ssize_t x, ssize_t y)
{
  // Carefully clone-and-modified over from float_image.c, but not
  // tested yet.
  //g_assert_not_reached ();

  // Reflect at image edges as advertised.
  if ( x < 0 ) {
    x = -x;
  }
  else if ( (size_t) x >= self->size_x ) {
    x = self->size_x - 2 - (x - self->size_x);
  }
  if ( y < 0 ) {
    y = -y;
  }
  else if ( (size_t) y >= self->size_y ) {
    y = self->size_y - 2 - (y - self->size_y);
  }

  return uint8_image_get_pixel (self, x, y);
}

void
uint8_image_statistics (UInt8Image *self, uint8_t *min, uint8_t *max,
                        double *mean, double *standard_deviation,
      gboolean use_mask_value, uint8_t mask_value)
{
  // Carefully clone-and-modified over from float_image.c, but not
  // tested yet.
  //g_assert_not_reached ();

  // Minimum and maximum sample values as integers.
  int imin = INT_MAX, imax = INT_MIN;

  // Buffer for one row of samples.
  uint8_t *row_buffer = g_new (uint8_t, self->size_x);

  *mean = 0.0;
  double s = 0.0;

  size_t sample_count = 0;      // Samples considered so far.
  size_t ii, jj;
  // If there is a mask value we are supposed to ignore,
  if ( use_mask_value ) {
    // iterate over all pixels, skipping pixels equal to mask value.
    for ( ii = 0 ; ii < self->size_y ; ii++ ) {
      asfPercentMeter((double)ii/(double)(self->size_y));
      uint8_image_get_row (self, ii, row_buffer);
      for ( jj = 0 ; jj < self->size_x ; jj++ ) {
        uint8_t cs = row_buffer[jj];   // Current sample.
        if ( cs == mask_value ) {
          continue;
        }
        if ( G_UNLIKELY (cs < imin) ) { imin = cs; }
        if ( G_UNLIKELY (cs > imax) ) { imax = cs; }
        double old_mean = *mean;
        *mean += (cs - *mean) / (sample_count + 1);
        s += (cs - old_mean) * (cs - *mean);
        sample_count++;
      }
    }
    asfPercentMeter(1.0);
  }
  else {
    // There is no mask value to ignore, so we do the same as the
    // above loop, but without the possible continue statement.
    for ( ii = 0 ; ii < self->size_y ; ii++ ) {
      asfPercentMeter((double)ii/(double)(self->size_y));
      uint8_image_get_row (self, ii, row_buffer);
      for ( jj = 0 ; jj < self->size_x ; jj++ ) {
        uint8_t cs = row_buffer[jj];   // Current sample.
        if ( G_UNLIKELY (cs < imin) ) { imin = cs; }
        if ( G_UNLIKELY (cs > imax) ) { imax = cs; }
        double old_mean = *mean;
        *mean += (cs - *mean) / (sample_count + 1);
        s += (cs - old_mean) * (cs - *mean);
        sample_count++;
      }
    }
    asfPercentMeter(1.0);
  }

  g_free (row_buffer);

  // Verify the new extrema have been found.
  g_assert (imin != INT_MAX);
  g_assert (imax != INT_MIN);

  // The new extrema had better be in the range supported by uint8_t.
  g_assert (imin >= 0);
  g_assert (imax <= UINT8_MAX);

  *min = imin;
  *max = imax;
  *standard_deviation = sqrt (s / (sample_count - 1));
}

int
uint8_image_band_statistics (UInt8Image *self, meta_stats *stats,
                             int line_count, int band_no,
                             gboolean use_mask_value, uint8_t mask_value)
{
  // Carefully cloned-and-modified over from float_image.c, but not
  // tested yet.
  //g_assert_not_reached ();

  // Minimum and maximum sample values as integers.
  int imin = INT_MAX, imax = INT_MIN;

  // Buffer for one row of samples.
  uint8_t *row_buffer = g_new (uint8_t, self->size_x);

  stats->mean = 0.0;
  double s = 0.0;

  size_t sample_count = 0;      // Samples considered so far.
  size_t ii, jj;
  // If there is a mask value we are supposed to ignore,
  if ( use_mask_value ) {
    // iterate over all pixels, skipping pixels equal to mask value.
    for ( ii = (band_no * line_count); // 0-ordered band number times lines is offset into image
          ii < (band_no+1) * line_count && ii < self->size_y;
          ii++ )
    {
      asfPercentMeter( (double)(ii - band_no*line_count)/(double)line_count );
      uint8_image_get_row (self, ii, row_buffer);
      for ( jj = 0 ; jj < self->size_x ; jj++ ) {
        uint8_t cs = row_buffer[jj];   // Current sample.
        if ( cs == mask_value ) {
          continue;
        }
        if ( G_UNLIKELY (cs < imin) ) { imin = cs; }
        if ( G_UNLIKELY (cs > imax) ) { imax = cs; }
        double old_mean = stats->mean;
        stats->mean += (cs - stats->mean) / (sample_count + 1);
        s += (cs - old_mean) * (cs - stats->mean);
        sample_count++;
      }
    }
    asfPercentMeter(1.0);
  }
  else {
    // There is no mask value to ignore, so we do the same as the
    // above loop, but without the possible continue statement.
    for ( ii = (band_no * line_count); // 0-ordered band number times lines is offset into image
          ii < (band_no+1) * line_count && ii < self->size_y;
          ii++ )
    {
      asfPercentMeter( (double)(ii - band_no*line_count)/(double)line_count );
      uint8_image_get_row (self, ii, row_buffer);
      for ( jj = 0 ; jj < self->size_x ; jj++ ) {
        uint8_t cs = row_buffer[jj];   // Current sample.
        if ( G_UNLIKELY (cs < imin) ) { imin = cs; }
        if ( G_UNLIKELY (cs > imax) ) { imax = cs; }
        double old_mean = stats->mean;
        stats->mean += (cs - stats->mean) / (sample_count + 1);
        s += (cs - old_mean) * (cs - stats->mean);
        sample_count++;
      }
    }
    asfPercentMeter(1.0);
  }

  g_free (row_buffer);

  // Verify the new extrema have been found.
  if (imin == INT_MAX || imax == INT_MIN)
    return 1;

  // The new extrema had better be in the range supported by uint8_t.
  if (imin < 0 || imax > UINT8_MAX)
    return 1;

  stats->min = imin;
  stats->max = imax;
  stats->std_deviation = sqrt (s / (sample_count - 1));

  return 0;
}

void
uint8_image_statistics_with_mask_interval (UInt8Image *self, uint8_t *min,
             uint8_t *max, double *mean,
             double *standard_deviation,
             uint8_t interval_start,
             uint8_t interval_end)
{
  // This method is a trivial clone-and-modify of
  // float_image_statistics, but it is totally untested at the moment.
  g_assert_not_reached ();

  // Minimum and maximum sample values as integers.
  int imin = INT_MAX, imax = INT_MIN;

  // Buffer for one row of samples.
  uint8_t *row_buffer = g_new (uint8_t, self->size_x);

  *mean = 0.0;
  double s = 0.0;

  size_t sample_count = 0;      // Samples considered so far.
  size_t ii, jj;
  for ( ii = 0 ; ii < self->size_y ; ii++ ) {
    uint8_image_get_row (self, ii, row_buffer);
    for ( jj = 0 ; jj < self->size_x ; jj++ ) {
      uint8_t cs = row_buffer[jj];   // Current sample.
      // If in the mask interval, do not consider this pixel any
      // further.
      if ( cs >= interval_start && cs <= interval_end ) {
  continue;
      }
      if ( G_UNLIKELY (cs < imin) ) { imin = cs; }
      if ( G_UNLIKELY (cs > imax) ) { imax = cs; }
      double old_mean = *mean;
      *mean += (cs - *mean) / (sample_count + 1);
      s += (cs - old_mean) * (cs - *mean);
      sample_count++;
    }
  }

  g_free (row_buffer);

  // Verify the new extrema have been found.
  g_assert (imin != INT_MAX);
  g_assert (imax != INT_MIN);

  // The new extrema had better be in the range supported by uint8_t.
  g_assert (imin >= 0);
  g_assert (imax <= UINT8_MAX);

  *min = imin;
  *max = imax;
  *standard_deviation = sqrt (s / (sample_count - 1));
}

void
uint8_image_approximate_statistics (UInt8Image *self, size_t stride,
                                    double *mean, double *standard_deviation,
            gboolean use_mask_value,
            uint8_t mask_value)
{
  // Rows and columns of samples that fit in image given stride
  // stride.
  size_t sample_columns = ceil (self->size_x / stride);
  size_t sample_rows = ceil (self->size_y / stride);
  // Total number of samples.
  size_t sample_count = sample_columns * sample_rows;

  // Create an image holding the sample values.
  UInt8Image *sample_image = uint8_image_new (sample_columns, sample_rows);

  // Load the sample values.
  size_t current_sample = 0;
  size_t ii;
  for ( ii = 0 ; ii < sample_columns ; ii++ ) {
    size_t jj;
    for ( jj = 0 ; jj < sample_rows ; jj++ ) {
      uint8_t sample = uint8_image_get_pixel (self, ii * stride, jj * stride);
      uint8_image_set_pixel (sample_image, ii, jj, sample);
      current_sample++;
    }
  }

  // Ensure that we got the right number of samples in our image.
  g_assert (current_sample == sample_count);

  // Compute the exact statistics of the sampled version of the image.
  // The _statistics method wants to compute min and max, so we let
  // it, even though we don't do anything with them (since they are
  // inaccurate).
  uint8_t min, max;
  uint8_image_statistics (sample_image, &min, &max, mean, standard_deviation,
                          use_mask_value, mask_value);

  uint8_image_free (sample_image);
}

void
uint8_image_approximate_statistics_with_mask_interval
  (UInt8Image *self, size_t stride, double *mean, double *standard_deviation,
   uint8_t interval_start, uint8_t interval_end)
{
  // This method is a trivial clone-and-modify of
  // float_image_approximate_statistics, but it is totally untested at
  // the moment.
  g_assert_not_reached ();

  // Rows and columns of samples that fit in image given stride
  // stride.
  size_t sample_columns = ceil (self->size_x / stride);
  size_t sample_rows = ceil (self->size_y / stride);
  // Total number of samples.
  size_t sample_count = sample_columns * sample_rows;

  // Create an image holding the sample values.
  UInt8Image *sample_image = uint8_image_new (sample_columns, sample_rows);

  // Load the sample values.
  size_t current_sample = 0;
  size_t ii;
  for ( ii = 0 ; ii < sample_columns ; ii++ ) {
    size_t jj;
    for ( jj = 0 ; jj < sample_rows ; jj++ ) {
      uint8_t sample = uint8_image_get_pixel (self, ii * stride, jj * stride);
      uint8_image_set_pixel (sample_image, ii, jj, sample);
      current_sample++;
    }
  }

  // Ensure that we got the right number of samples in our image.
  g_assert (current_sample == sample_count);

  // Compute the exact statistics of the sampled version of the image.
  // The _statistics method wants to compute min and max, so we let
  // it, even though we don't do anything with them (since they are
  // inaccurate).
  uint8_t min, max;
  uint8_image_statistics_with_mask_interval (sample_image, &min, &max,
               mean, standard_deviation,
               interval_start, interval_end);

  uint8_image_free (sample_image);
}

gsl_histogram *
uint8_image_gsl_histogram (UInt8Image *self, double min, double max,
                           size_t bin_count)
{
  // Carefully clone-and-modified over from float_image.c, but not
  // tested yet.
  g_assert_not_reached ();

  // Initialize the histogram.
  gsl_histogram *histogram = gsl_histogram_alloc (bin_count);
  gsl_histogram_set_ranges_uniform (histogram, min, max);

  // Buffer for one row of samples.
  uint8_t *row_buffer = g_new (uint8_t, self->size_x);

  // Populate the histogram over every sample in the image.
  size_t ii, jj;
  for (ii = 0 ; ii < self->size_y ; ii++ ) {
    uint8_image_get_row (self, ii, row_buffer);
    for ( jj = 0 ; jj < self->size_x ; jj++ ) {
      gsl_histogram_increment (histogram, row_buffer[jj]);
    }
  }

  g_free (row_buffer);

  return histogram;
}

double
uint8_image_apply_kernel (UInt8Image *self, ssize_t x, ssize_t y,
        gsl_matrix *kern)
{
  // Carefully clone-and-modified over from float_image.c, but not
  // tested yet.
  g_assert_not_reached ();

  g_assert (x >= 0 && (size_t) x < self->size_x);
  g_assert (y >= 0 && (size_t) y < self->size_y);
  g_assert (kern->size2 % 2 == 1);
  g_assert (kern->size2 == kern->size1);

  size_t ks = kern->size2;    // Kernel size.

  double sum = 0;                // Result.

  size_t ii;
  for ( ii = 0 ; ii < kern->size1 ; ii++ ) {
    ssize_t iy = y - ks / 2 + ii; // Current image y pixel index.
    size_t jj;
    for ( jj = 0 ; jj < kern->size2 ; jj++ ) {
      ssize_t ix = x - ks / 2 + jj; // Current image x pixel index
      sum += (gsl_matrix_get (kern, jj, ii)
              * uint8_image_get_pixel_with_reflection (self, ix, iy));
    }
  }

  return sum;
}

double
uint8_image_sample (UInt8Image *self, double x, double y,
                    uint8_image_sample_method_t sample_method)
{
  g_assert (x >= 0.0 && x <= (double) self->size_x - 1.0);
  g_assert (y >= 0.0 && y <= (double) self->size_y - 1.0);

  switch ( sample_method ) {

  case UINT8_IMAGE_SAMPLE_METHOD_NEAREST_NEIGHBOR:
    return uint8_image_get_pixel (self, round (x), round (y));
    break;

  case UINT8_IMAGE_SAMPLE_METHOD_BILINEAR:
    {
      // Indicies of points we are interpolating between (x below, y
      // below, etc., where below is interpreted in the numerical
      // sense, not the image orientation sense.).
      size_t xb = floor (x), yb = floor (y), xa = ceil (x), ya = ceil (y);
      size_t ts = self->tile_size;   // Convenience alias.
      // Offset of xb, yb, etc. relative to tiles they lie in.
      size_t xbto = xb % ts, ybto = yb % ts, xato = xa % ts, yato = ya % ts;
      // Values of points we are interpolating between.
      uint8_t ul, ur, ll, lr;

      // If the points were are interpolating between don't span a
      // tile edge, we load them straight from tile memory to save
      // some time.
      if ( G_LIKELY (   xbto != ts - 1 && xato != 0
                     && ybto != ts - 1 && yato != 0) ) {
        // The tile indicies.
        size_t tx = xb / ts, ty = yb / ts;
        // Tile offset in flattened list of tile addresses.
        size_t tile_offset = ty * self->tile_count_x + tx;
        uint8_t *tile_address = self->tile_addresses[tile_offset];
        if ( G_UNLIKELY (tile_address == NULL) ) {
          tile_address = load_tile (self, tx, ty);
        }
        ul = tile_address[ybto * self->tile_size + xbto];
        ur = tile_address[ybto * self->tile_size + xato];
        ll = tile_address[yato * self->tile_size + xbto];
        lr = tile_address[yato * self->tile_size + xato];
      }
      else {
        // We are spanning a tile edge, so we just get the pixels
        // using the inefficient but easy get_pixel method.
        ul = uint8_image_get_pixel (self, floor (x), floor (y));
        ur = uint8_image_get_pixel (self, ceil (x), floor (y));
        ll = uint8_image_get_pixel (self, floor (x), ceil (y));
        lr = uint8_image_get_pixel (self, ceil (x), ceil (y));
      }

      // Upper and lower values interpolated in the x direction.
      double ux = ul + (ur - ul) * (x - floor (x));
      double lx = ll + (lr - ll) * (x - floor (x));

      return ux + (lx - ux) * (y - floor (y));
    }
    break;
  case UINT8_IMAGE_SAMPLE_METHOD_BICUBIC:
    {
      // Should never be here ...bicubic resampling can result in negative
      // values and should not be used for resampling unsigned values...
      //g_assert_not_reached ();
      asfPrintError ("BICUBIC resampling for BYTE data is not supported.\n");

      static gboolean first_time_through = TRUE;
      // Splines in the x direction, and their lookup accelerators.
      static double *x_indicies;
      static double *values;
      static gsl_spline **xss;
      static gsl_interp_accel **xias;
      // Spline between splines in the y direction, and lookup accelerator.
      static double *y_spline_indicies;
      static double *y_spline_values;
      static gsl_spline *ys;
      static gsl_interp_accel *yia;

      // All these splines have size 4.
      const size_t ss = 4;

      size_t ii;                // Index variable.

      if ( first_time_through ) {
        // Allocate memory for the splines in the x direction.
        x_indicies = g_new (double, ss);
        values = g_new (double, ss);
        xss = g_new (gsl_spline *, ss);
        xias = g_new (gsl_interp_accel *, ss);
        for ( ii = 0 ; ii < ss ; ii++ ) {
          xss[ii] = gsl_spline_alloc (gsl_interp_cspline, ss);
          xias[ii] = gsl_interp_accel_alloc ();
        }

        // Allocate memory for the spline in the y direction.
        y_spline_indicies = g_new (double, ss);
        y_spline_values = g_new (double, ss);
        ys = gsl_spline_alloc (gsl_interp_cspline, ss);
        yia = gsl_interp_accel_alloc ();
        first_time_through = FALSE;
      }

      // Get the values for the nearest 16 points.
      size_t jj;                // Index variable.
      for ( ii = 0 ; ii < ss ; ii++ ) {
        for ( jj = 0 ; jj < ss ; jj++ ) {
          x_indicies[jj] = floor (x) - 1 + jj;
          values[jj]
            = uint8_image_get_pixel_with_reflection (self, x_indicies[jj],
                                                     floor (y) - 1 + ii);
        }
        gsl_spline_init (xss[ii], x_indicies, values, ss);
      }

      // Set up the spline that runs in the y direction.
      for ( ii = 0 ; ii < ss ; ii++ ) {
        y_spline_indicies[ii] = floor (y) - 1 + ii;
        y_spline_values[ii] = gsl_spline_eval (xss[ii], x, xias[ii]);
      }
      gsl_spline_init (ys, y_spline_indicies, y_spline_values, ss);

      double ret_val = gsl_spline_eval (ys, y, yia);
// NOTE... NOTE... BICUBIC resample returns negative values if the byte values are
// too close to zero and the spline fit goes negative in the neighborhood of the pixel
/*
      if (ret_val > 255.0 || ret_val < 0.0) {
        asfPrintWarning("Bicubic resampling of BYTE data returned out of range value (%f).\n"
          "...Continuing, but negative values will be forced to zero, and values above 255\n"
          " will be forced to 255\n",
          ret_val);
        if (ret_val > 265.0) {
          asfPrintError("Bicubic resampling of BYTE data returned a value (%f) too far above 255.0 to\n"
            "cap to 255.0\n", ret_val);
        }
        if (ret_val < -10.0) {
          asfPrintError("Bicubic resampling of BYTE data returned a value (%f) too far below 0.0 to\n"
            "cap to 0.0\n", ret_val);
        }
        ret_val = ret_val < 0.0 ? 0.0 : ret_val;
        ret_val = ret_val > 255.0 ? 255.0 : ret_val;
      }
*/
      return ret_val;
    }
    break;
  default:
    g_assert_not_reached ();
    return -42;         // Reassure the compiler.
  }
}

gboolean
uint8_image_equals (UInt8Image *self, UInt8Image *other)
{
  // Compare image sizes.
  if ( self->size_x != other->size_x ) {
    return FALSE;
  }
  if ( self->size_y != other->size_y ) {
    return FALSE;
  }

  size_t sz = self->size_x; // Convenience alias.

  // Compare image pixels.
  size_t ii, jj;
  for ( ii = 0 ; ii < sz ; ii++ ) {
    for ( jj = 0 ; jj < sz ; jj++ ) {
      if ( G_UNLIKELY (uint8_image_get_pixel (self, jj, ii)
           != uint8_image_get_pixel (other, jj, ii)) ) {
  return FALSE;
      }
    }
  }

  return TRUE;
}

// Flip an image about a horizontal line through the center of the image
void
uint8_image_flip_y(UInt8Image *self)
{
  size_t ii, jj;

  asfLineMeter(jj, self->size_y);
  for (jj = 0; jj < self->size_y / 2; ++jj) {
    asfLineMeter(2 * jj + 1, self->size_y);
    size_t jj2 = self->size_y - 1 - jj;
    for (ii = 0; ii < self->size_x; ++ii) {
      uint8 a = uint8_image_get_pixel(self, ii, jj);
      uint8 b = uint8_image_get_pixel(self, ii, jj2);
      uint8_image_set_pixel(self, ii, jj, b);
      uint8_image_set_pixel(self, ii, jj2, a);
    }
  }
  asfLineMeter(1, 1);
}

// Flip an image about a vertical line through the center of the image
void
uint8_image_flip_x(UInt8Image *self)
{
  size_t ii, jj;

  for (ii = 0; ii < self->size_x / 2; ++ii) {
    asfLineMeter(2 * ii + 1, self->size_y);
    size_t ii2 = self->size_x - 1 - ii;
    for (jj = 0; jj < self->size_y; ++jj) {
      uint8 a = uint8_image_get_pixel(self, ii, jj);
      uint8 b = uint8_image_get_pixel(self, ii2, jj);
      uint8_image_set_pixel(self, ii, jj, b);
      uint8_image_set_pixel(self, ii2, jj, a);
    }
  }
  asfLineMeter(1, 1);
}

// Bring the tile cache file on the disk fully into sync with the
// latest image data stored in the memory cache.
static void
synchronize_tile_file_with_memory_cache (UInt8Image *self)
{
  // If we aren't using a tile file, this operation doesn't make
  // sense.
  g_assert (self->tile_file != NULL);

  guint ii;
  for ( ii = 0 ; ii < self->tile_queue->length ; ii++ ) {
    size_t tile_offset = GPOINTER_TO_INT (g_queue_peek_nth (self->tile_queue,
                  ii));
    cached_tile_to_disk (self, tile_offset);
  }
}

void
uint8_image_freeze (UInt8Image *self, FILE *file_pointer)
{
  FILE *fp = file_pointer;  // Convenience alias.

  g_assert (file_pointer != NULL);

  size_t write_count = fwrite (&(self->size_x), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->size_y), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->cache_space), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->cache_area), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->tile_size), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->cache_size_in_tiles), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->tile_count_x), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->tile_count_y), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->tile_count), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  write_count = fwrite (&(self->tile_area), sizeof (size_t), 1, fp);
  g_assert (write_count == 1);

  // We don't bother serializing the cache -- its a pain to keep track
  // of and probably almost never worth it.

  // We write the tile queue pointer away, so that when we later thaw
  // the serialized version, we can tell if a cache file is in use or
  // not (if it isn't tile_queue will be NULL).
  write_count = fwrite (&(self->tile_queue), sizeof (GQueue *), 1, fp);
  g_assert (write_count == 1);

  // If there was no cache file...
  if ( self->tile_queue == NULL ) {
    // We store the contents of the first tile and are done.
    write_count = fwrite (self->tile_addresses[0], sizeof (uint8_t),
        self->tile_area, fp);
    if ( write_count < self->tile_area ) {
      if ( ferror (fp) ) {
  fprintf (stderr, "Error writing serialized UInt8Image instance during "
     "freeze: %s\n", strerror (errno));
  exit (EXIT_FAILURE);
      }
    }
    g_assert (write_count == self->tile_area);
  }
  // otherwise, the in memory cache needs to be copied into the tile
  // file and the tile file saved in the serialized version of self.
  else {
    synchronize_tile_file_with_memory_cache (self);
    uint8_t *buffer = g_new (uint8_t, self->tile_area);
    size_t ii;
    off_t tmp = FTELL64 (self->tile_file);
    int return_code = FSEEK64 (self->tile_file, 0, SEEK_SET);
    g_assert (return_code == 0);
    for ( ii = 0 ; ii < self->tile_count ; ii++ ) {
      size_t read_count = fread (buffer, sizeof (uint8_t), self->tile_area,
         self->tile_file);
      g_assert (read_count == self->tile_area);
      write_count = fwrite (buffer, sizeof (uint8_t), self->tile_area, fp);
      g_assert (write_count == self->tile_area);
    }
    return_code = FSEEK64 (self->tile_file, tmp, SEEK_SET);
    g_assert (return_code == 0);
    g_free (buffer);
  }
}

int
uint8_image_band_store(UInt8Image *self, const char *file,
           meta_parameters *meta, int append_flag)
{
  // Give status
  if (meta->general->band_count == 1)
    asfPrintStatus("\n\nStoring image ...\n");
  else
    asfPrintStatus("\n\nStoring band ...\n");

  // Open the file to write to.
  FILE *fp = fopen (file, append_flag ? "ab" : "wb");
  // FIXME: we need some error handling and propagation here.
  g_assert (fp != NULL);

  // We will write the image data in horizontal stips one line at a time.
  uint8_t *line_buffer = g_new (uint8_t, self->size_x);

  // Sanity check
  if (meta->general->line_count != (int)self->size_y ||
      meta->general->sample_count != (int)self->size_x)
  {
      asfPrintError("Inconsistency between metadata and image!\n"
                    "Metadata says: %dx%d LxS, image has %dx%d\n"
                    "Possibly did not write metadata before storing image.\n",
                    meta->general->line_count, meta->general->sample_count,
                    self->size_y, self->size_x);
  }

  // Reorganize data into tiles in tile oriented disk file.
  int ii;
  for ( ii = 0 ; ii < (int)self->size_y ; ii++ ) {
    uint8_image_get_row (self, ii, line_buffer);

    size_t write_count =
      fwrite(line_buffer, sizeof(uint8_t), self->size_x, fp);

    if (write_count < self->size_x) {
      // it must have been a write error (no space left, possibly)
      g_assert(ferror(self->tile_file));
      // so print an error message
      fprintf(stderr, "Error writing file %s: %s\n", file, strerror(errno));
      // and exit
      exit (EXIT_FAILURE);
    }

    g_assert(write_count == self->size_x);
  }

  // Done with the line buffer.
  g_free (line_buffer);

  // Close file being written.
  int return_code = fclose (fp);
  g_assert (return_code == 0);

  // Return success code.
  return 0;
}

int
uint8_image_store (UInt8Image *self, const char *file)
{
  meta_parameters *meta;
  meta = meta_read(file);

  int ret = uint8_image_band_store(self, file, meta, 0);
  meta_free(meta);

  return ret;
}

int
uint8_image_export_as_jpeg (UInt8Image *self, const char *file,
                            size_t max_dimension, gboolean use_mask_value,
          uint8_t mask_value)
{
  // Carefully clone-and-modified over from float_image.c, but not
  // tested yet.
  g_assert_not_reached ();

  size_t scale_factor;          // Scale factor to use for output image.
  if ( self->size_x > self->size_y ) {
    scale_factor = ceil ((double) self->size_x / max_dimension);
  }
  else {
    scale_factor = ceil ((double) self->size_y / max_dimension);
  }

  // We want the scale factor to be odd, so that we can easily use a
  // standard kernel to average things.
  if ( scale_factor % 2 == 0 ) {
    scale_factor++;
  }

  // Output JPEG x and y dimensions.
  size_t osx = self->size_x / scale_factor;
  size_t osy = self->size_y / scale_factor;

  // Number of pixels in output image.
  size_t pixel_count = osx * osy;

  // Pixels of the output image.
  unsigned char *pixels = g_new (unsigned char, pixel_count);

  JSAMPLE test_jsample;         // For verifying properties of JSAMPLE type.
  /* Here are some very funky checks to try to ensure that the JSAMPLE
     really is the type we expect, so we can scale properly.  */
  g_assert (sizeof (unsigned char) == 1);
  g_assert (sizeof (unsigned char) == sizeof (JSAMPLE));
  test_jsample = 0;
  test_jsample--;
  g_assert (test_jsample == UCHAR_MAX);

  // Stuff needed by libjpeg.
  struct jpeg_compress_struct cinfo;
  struct jpeg_error_mgr jerr;
  cinfo.err = jpeg_std_error (&jerr);
  jpeg_create_compress (&cinfo);

  // Open output file.
  FILE *fp = fopen (file, "wb");
  if ( fp == NULL ) { perror ("error opening file"); }
  // FIXME: we need some error handling and propagation here.
  g_assert (fp != NULL);

  // Connect jpeg output to the output file to be used.
  jpeg_stdio_dest (&cinfo, fp);

  // Set image parameters that libjpeg needs to know about.
  cinfo.image_width = osx;
  cinfo.image_height = osy;
  cinfo.input_components = 1;   // Grey scale => 1 color component / pixel.
  cinfo.in_color_space = JCS_GRAYSCALE;
  jpeg_set_defaults (&cinfo);   // Use default compression parameters.
  // Reassure libjpeg that we will be writing a complete JPEG file.
  jpeg_start_compress (&cinfo, TRUE);

  // As advertised, we will average pixels together.
  g_assert (scale_factor % 2 != 0);
  size_t kernel_size = scale_factor;
  gsl_matrix *averaging_kernel
    = gsl_matrix_alloc (kernel_size, kernel_size);
  double kernel_value = 1.0 / ((double)kernel_size * kernel_size);
  size_t ii, jj;                // Index values.
  for ( ii = 0 ; ii < averaging_kernel->size1 ; ii++ ) {
    for ( jj = 0 ; jj < averaging_kernel->size2 ; jj++ ) {
      gsl_matrix_set (averaging_kernel, ii, jj, kernel_value);
    }
  }

  // Sample input image, putting scaled results into output image.
  size_t sample_stride = scale_factor;
  for ( ii = 0 ; ii < osy ; ii++ ) {
    for ( jj = 0 ; jj < osx ; jj++ ) {
      // Input image average pixel value.
      double ival = uint8_image_apply_kernel (self, jj * sample_stride,
                ii * sample_stride,
                averaging_kernel);
      // Set output value.
      int32_t oval = round (ival);
      // In case floating point arithmetic wierdness gets us in
      // trouble, we correct.
      if ( oval < 0 ) {
  oval = 0;
      }
      else if ( oval > UINT8_MAX ) {
  oval = UINT8_MAX;
      }
      pixels[ii * osx + jj] = oval;
    }
  }

  // Write the jpeg, one row at a time.
  const int rows_to_write = 1;
  JSAMPROW *row_pointer = g_new (JSAMPROW, rows_to_write);
  while ( cinfo.next_scanline < cinfo.image_height ) {
    int rows_written;
    row_pointer[0] = &(pixels[cinfo.next_scanline * osx]);
    rows_written = jpeg_write_scanlines (&cinfo, row_pointer, rows_to_write);
    g_assert (rows_written == rows_to_write);
  }
  g_free (row_pointer);

  // Finsh compression and close the jpeg.
  jpeg_finish_compress (&cinfo);
  int return_code = fclose (fp);
  g_assert (return_code == 0);
  jpeg_destroy_compress (&cinfo);

  g_free (pixels);

  return 0;                     // Return success indicator.
}

size_t
uint8_image_get_cache_size (UInt8Image *self)
{
  g_assert_not_reached ();      // Stubbed out for now.
  // Compiler reassurance.
  self = self;
  return 0;
}

void
uint8_image_set_cache_size (UInt8Image *self, size_t size)
{
  g_assert_not_reached ();      // Stubbed out for now.
  // Compiler reassurance.
  self = self; size = size;
}

void
uint8_image_free (UInt8Image *self)
{
  // Close the tile file (which shouldn't have to remove it since its
  // already unlinked), if we were ever using it.
  if ( self->tile_file != NULL ) {
    int return_code = fclose (self->tile_file);
    g_assert (return_code == 0);
  }

  // Deallocate dynamic memory.

  g_free (self->tile_addresses);

  // If we didn't need a tile file, we also won't have a tile queue.
  if ( self->tile_queue != NULL ) {
    g_queue_free (self->tile_queue);
  }

  g_free (self->cache);

  if (self->tile_file_name) {
#ifdef win32
     unlink_tmp_file(self->tile_file_name->str);
     g_string_free(self->tile_file_name, TRUE);
#endif
  }

  g_free (self);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#pragma once

#include <memory>

#include <gsl/span>

#include "graphics/math.h"
#include "infrastructure/macros.h"

namespace gfx {
	using MdfRenderMaterialPtr = std::shared_ptr<class MdfRenderMaterial>;

	class RenderingDevice;
	class Texture;
	using TextureRef = std::shared_ptr<Texture>;

	struct Line2d {
		XMFLOAT2 from;
		XMFLOAT2 to;
		XMCOLOR diffuse;

		Line2d() {}

		Line2d(XMFLOAT2 from, XMFLOAT2 to, XMCOLOR color)
			: from(from), to(to), diffuse(color) {}
	};

#pragma pack(push, 1)
	struct Vertex2d {
		XMFLOAT4 pos;
		XMFLOAT4 normal;
		XMCOLOR diffuse;
		XMFLOAT2 uv;
	};
#pragma pack(pop)

	enum class SamplerType2d {
		CLAMP,
		WRAP,
		POINT
	};

	class SamplerState;

	/*
		Renders shapes in 2d screen space.
	*/
	class ShapeRenderer2d {
	public:
		explicit ShapeRenderer2d(RenderingDevice& g);
		~ShapeRenderer2d();

		void DrawRectangle(float x, float y, float width, float height, gfx::Texture& texture, uint32_t color = 0xFFFFFFFF, SamplerType2d samplerType = SamplerType2d::CLAMP) {
			DrawRectangle(x, y, width, height, &texture, color, samplerType);
		}

		void DrawRectangle(float x, float y, float width, float height, uint32_t color) {
			DrawRectangle(x, y, width, height, nullptr, color);
		}

		void DrawRectangle(gsl::span<Vertex2d, 4> corners,
			gfx::Texture* texture,
			gfx::Texture* mask = nullptr,
			SamplerType2d samplerType = SamplerType2d::CLAMP,
			bool blending = true);

		void DrawRectangle(gsl::span<Vertex2d, 4> corners,
			const gfx::MdfRenderMaterialPtr &material);

		void DrawRectangle(gsl::span<Vertex2d, 4> corners);

		void DrawLines(gsl::span<Line2d> lines);

		void DrawRectangleOutline(XMFLOAT2 topLeft, XMFLOAT2 bottomRight, XMCOLOR color);

		void DrawFullScreenQuad();

		/**
		 * Renders a circle/pie segment for use with the radial menu.
  		 */
		void DrawPieSegment(int segments,
			int x, int y,
			float angleCenter, float angleWidth,
			int innerRadius, int innerOffset,
			int outerRadius, int outerOffset,
			XMCOLOR color1, XMCOLOR color2);

		NO_COPY_OR_MOVE(ShapeRenderer2d);

	private:

		void DrawRectangle(
			float x, float y, float width, float height,
			gfx::Texture* texture,
			uint32_t color = 0xFFFFFFFF,
			SamplerType2d samplerType = SamplerType2d::CLAMP
		);

		SamplerState &getSamplerState(SamplerType2d type) const;

		struct Impl;
		std::unique_ptr<Impl> mImpl;
	};

}

ext: h

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
#ifndef _UTIL_H
#define _UTIL_H

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <chrono>
#include <string>
#include <random>
#include <algorithm>
#include <mutex>
#include <atomic>
#include <tbb/pipeline.h>
#include "blocks.pb.h"
#ifdef __APPLE__
extern "C"
{
#include <cblas.h>
}
extern "C"
{
  void cblas_saxpy(const int N, const float alpha, const float *X,
		   const int incX, float *Y, const int incY);
  void cblas_scopy(const int N, const float *X, const int incX,
		   float *Y, const int incY);
  float cblas_sdot(const int N, const float  *X, const int incX,
		    const float  *Y, const int incY);
}
#else
#include "mkl.h"
#endif
typedef unsigned long long uint64;
typedef unsigned int       uint32;
typedef std::chrono::high_resolution_clock Time;
#ifdef LEVEL1_DCACHE_LINESIZE
#define CACHE_LINE_SIZE LEVEL1_DCACHE_LINESIZE
#else
#define CACHE_LINE_SIZE 64
#endif

typedef struct {
    int u_, v_;
    float r_;
} Record;

extern std::chrono::time_point<Time> s,e;
extern std::default_random_engine generator;
extern std::normal_distribution<float> gaussian;

#ifdef FETCH
inline void prefetch_range(char *addr, size_t len) {
  char *cp;
  char *end = addr + len;

  for (cp = addr; cp < end; cp += CACHE_LINE_SIZE)
    __builtin_prefetch(cp,1,0);
}
#endif

inline void align_alloc(float** u, int nu, int dim) {
    int piece = nu/1050000+1;
    int nn = nu/piece;
    int k;
    for(k=0; k<piece-1; k++) {
        u[k*nn] = (float*)mkl_malloc(nn*dim*sizeof(float), CACHE_LINE_SIZE);
        for(int i=1; i<nn; i++)
            u[k*nn+i] = u[k*nn+i-1] + dim;
    }
    u[k*nn] = (float*)mkl_malloc((nn+nu%piece)*dim*sizeof(float), CACHE_LINE_SIZE);
    for(int i=1; i<nn+nu%piece; i++)
        u[k*nn+i] = u[k*nn+i-1] + dim;
}

inline void plain_read(const char* data, mf::Blocks& blocks) {
  FILE* fr = fopen(data, "rb");
  std::vector<char> buf;
  uint32 isize;
  mf::Block* bk;
  while(fread(&isize, 1, sizeof(isize), fr)) {
      buf.resize(isize);
      fread((char*)buf.data(), 1, isize, fr);
      bk = blocks.add_block();
      bk->ParseFromArray(buf.data(), isize);
  }
  fclose(fr);
}

inline float active(float val, int type) {
  switch(type) {
    case 0: return val;                     //least square
    case 1: return 1.0f/(1.0f+expf(-val));  //sigmoid
  }
}
inline float cal_grad(float r, float pred, int type) {
  switch(type) {
    case 0: return r - pred;          //least square
    case 1: return r - pred;          //0-1 logistic regression
  }
}

inline float next_float(){
    return static_cast<float>( rand() ) / (static_cast<float>( RAND_MAX )+1.0);
}

inline float next_float2(){
    return (static_cast<float>( rand() ) + 1.0 ) / (static_cast<float>(RAND_MAX) + 2.0);
}

inline float normsqr(float* x, int num) {
   return cblas_sdot(num, x, 1, x, 1);
}

inline float sample_normal(){
    float x,y,s;
    do{
        x = 2 * next_float2() - 1.0;
        y = 2 * next_float2() - 1.0;
        s = x*x + y*y;
    }while( s >= 1.0 || s == 0.0 );

    return x * sqrt( -2.0 * log(s) / s ) ;
}

inline float sample_gamma( float alpha, float beta ) {
    if ( alpha < 1.0 ) {
        float u;
        do {
            u = next_float();
        } while (u == 0.0);
        return sample_gamma(alpha + 1.0, beta) * pow(u, 1.0 / alpha);
    } else {
        float d,c,x,v,u;
        d = alpha - 1.0/3.0;
        c = 1.0 / sqrt( 9.0 * d );
        do {
            do {
                x = sample_normal();
                v = 1.0 + c*x;
            } while ( v <= 0.0 );
            v = v * v * v;
            u = next_float();
        } while ( (u >= (1.0 - 0.0331 * (x*x) * (x*x)))
            && (log(u) >= (0.5 * x * x + d * (1.0 - v + log(v)))) );
        return d * v / beta;
    }
}

inline void gamma_posterior( float &lambda, float prior_alpha, float prior_beta, float psum_sqr, float psum_cnt ){
    float alpha = prior_alpha + 0.5*psum_cnt;
    float beta  = prior_beta + 0.5*psum_sqr;
    lambda = sample_gamma( alpha, beta );
}

inline void normsqr_col(float** m, int d, int size, float* norm) {
#pragma omp parallel for
    for(int i=0; i<d; i++) {
        for(int j=0; j<size; j++) norm[i] += m[j][i]*m[j][i];
    }
}

inline int padding(int dim) {
    return ((dim*sizeof(float)-1)/CACHE_LINE_SIZE*CACHE_LINE_SIZE+CACHE_LINE_SIZE)/sizeof(float);
}

#endif

max_stars_repo_name: dmlc/mf

max_issues_repo_name: dmlc/mf

max_forks_repo_name: dmlc/mf


text:
#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>

int
main (void)
{
  float f = 1.0/3.0;
  double d = 1.0/3.0;

  double fd = f; /* promote from float to double */

  printf (" f="); gsl_ieee_printf_float(&f);
  printf ("\n");

  printf ("fd="); gsl_ieee_printf_double(&fd);
  printf ("\n");

  printf (" d="); gsl_ieee_printf_double(&d);
  printf ("\n");

  return 0;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef libceed_petsc_examples_setup_h
#define libceed_petsc_examples_setup_h

#include <ceed.h>
#include <petsc.h>

#include "structs.h"

PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data);
PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
                                    CeedInt topo_dim, CeedInt q_extra,
                                    PetscInt num_comp_x, PetscInt num_comp_u,
                                    PetscInt g_size, PetscInt xl_size,
                                    BPData bp_data, CeedData data,
                                    PetscBool setup_rhs, CeedVector rhs_ceed,
                                    CeedVector *target);
PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt num_levels,
                                      CeedInt num_comp_u, CeedData *data, CeedInt *leveldegrees,
                                      CeedQFunction qf_restrict, CeedQFunction qf_prolong);

#endif // libceed_petsc_examples_setup_h

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:44:49 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -hc2hc-forward 10 */

/*
 * This function contains 168 FP additions, 84 FP multiplications,
 * (or, 126 additions, 42 multiplications, 42 fused multiply/add),
 * 43 stack variables, and 80 memory accesses
 */
static const fftw_real K587785252 = FFTW_KONST(+0.587785252292473129168705954639072768597652438);
static const fftw_real K951056516 = FFTW_KONST(+0.951056516295153572116439333379382143405698634);
static const fftw_real K250000000 = FFTW_KONST(+0.250000000000000000000000000000000000000000000);
static const fftw_real K559016994 = FFTW_KONST(+0.559016994374947424102293417182819058860154590);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_hc2hc_forward_10(fftw_real *A, const fftw_complex *W, int iostride, int m, int dist)
{
     int i;
     fftw_real *X;
     fftw_real *Y;
     X = A;
     Y = A + (10 * iostride);
     {
	  fftw_real tmp170;
	  fftw_real tmp181;
	  fftw_real tmp162;
	  fftw_real tmp175;
	  fftw_real tmp165;
	  fftw_real tmp176;
	  fftw_real tmp166;
	  fftw_real tmp183;
	  fftw_real tmp155;
	  fftw_real tmp178;
	  fftw_real tmp158;
	  fftw_real tmp179;
	  fftw_real tmp159;
	  fftw_real tmp182;
	  fftw_real tmp168;
	  fftw_real tmp169;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp168 = X[0];
	  tmp169 = X[5 * iostride];
	  tmp170 = tmp168 - tmp169;
	  tmp181 = tmp168 + tmp169;
	  {
	       fftw_real tmp160;
	       fftw_real tmp161;
	       fftw_real tmp163;
	       fftw_real tmp164;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp160 = X[4 * iostride];
	       tmp161 = X[9 * iostride];
	       tmp162 = tmp160 - tmp161;
	       tmp175 = tmp160 + tmp161;
	       tmp163 = X[6 * iostride];
	       tmp164 = X[iostride];
	       tmp165 = tmp163 - tmp164;
	       tmp176 = tmp163 + tmp164;
	  }
	  tmp166 = tmp162 + tmp165;
	  tmp183 = tmp175 + tmp176;
	  {
	       fftw_real tmp153;
	       fftw_real tmp154;
	       fftw_real tmp156;
	       fftw_real tmp157;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp153 = X[2 * iostride];
	       tmp154 = X[7 * iostride];
	       tmp155 = tmp153 - tmp154;
	       tmp178 = tmp153 + tmp154;
	       tmp156 = X[8 * iostride];
	       tmp157 = X[3 * iostride];
	       tmp158 = tmp156 - tmp157;
	       tmp179 = tmp156 + tmp157;
	  }
	  tmp159 = tmp155 + tmp158;
	  tmp182 = tmp178 + tmp179;
	  {
	       fftw_real tmp167;
	       fftw_real tmp171;
	       fftw_real tmp172;
	       fftw_real tmp186;
	       fftw_real tmp184;
	       fftw_real tmp185;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp167 = K559016994 * (tmp159 - tmp166);
	       tmp171 = tmp159 + tmp166;
	       tmp172 = tmp170 - (K250000000 * tmp171);
	       X[iostride] = tmp167 + tmp172;
	       X[3 * iostride] = tmp172 - tmp167;
	       X[5 * iostride] = tmp170 + tmp171;
	       tmp186 = K559016994 * (tmp182 - tmp183);
	       tmp184 = tmp182 + tmp183;
	       tmp185 = tmp181 - (K250000000 * tmp184);
	       X[2 * iostride] = tmp185 - tmp186;
	       X[4 * iostride] = tmp186 + tmp185;
	       X[0] = tmp181 + tmp184;
	  }
	  {
	       fftw_real tmp173;
	       fftw_real tmp174;
	       fftw_real tmp177;
	       fftw_real tmp180;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp173 = tmp155 - tmp158;
	       tmp174 = tmp162 - tmp165;
	       Y[-iostride] = -((K951056516 * tmp173) + (K587785252 * tmp174));
	       Y[-3 * iostride] = (K587785252 * tmp173) - (K951056516 * tmp174);
	       tmp177 = tmp175 - tmp176;
	       tmp180 = tmp178 - tmp179;
	       Y[-2 * iostride] = (K951056516 * tmp177) - (K587785252 * tmp180);
	       Y[-4 * iostride] = (K951056516 * tmp180) + (K587785252 * tmp177);
	  }
     }
     X = X + dist;
     Y = Y - dist;
     for (i = 2; i < m; i = i + 2, X = X + dist, Y = Y - dist, W = W + 9) {
	  fftw_real tmp39;
	  fftw_real tmp87;
	  fftw_real tmp132;
	  fftw_real tmp144;
	  fftw_real tmp73;
	  fftw_real tmp84;
	  fftw_real tmp85;
	  fftw_real tmp91;
	  fftw_real tmp92;
	  fftw_real tmp93;
	  fftw_real tmp100;
	  fftw_real tmp103;
	  fftw_real tmp128;
	  fftw_real tmp121;
	  fftw_real tmp122;
	  fftw_real tmp142;
	  fftw_real tmp50;
	  fftw_real tmp61;
	  fftw_real tmp62;
	  fftw_real tmp88;
	  fftw_real tmp89;
	  fftw_real tmp90;
	  fftw_real tmp107;
	  fftw_real tmp110;
	  fftw_real tmp127;
	  fftw_real tmp118;
	  fftw_real tmp119;
	  fftw_real tmp141;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp33;
	       fftw_real tmp131;
	       fftw_real tmp38;
	       fftw_real tmp130;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp33 = X[0];
	       tmp131 = Y[-9 * iostride];
	       {
		    fftw_real tmp35;
		    fftw_real tmp37;
		    fftw_real tmp34;
		    fftw_real tmp36;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp35 = X[5 * iostride];
		    tmp37 = Y[-4 * iostride];
		    tmp34 = c_re(W[4]);
		    tmp36 = c_im(W[4]);
		    tmp38 = (tmp34 * tmp35) - (tmp36 * tmp37);
		    tmp130 = (tmp36 * tmp35) + (tmp34 * tmp37);
	       }
	       tmp39 = tmp33 - tmp38;
	       tmp87 = tmp33 + tmp38;
	       tmp132 = tmp130 + tmp131;
	       tmp144 = tmp131 - tmp130;
	  }
	  {
	       fftw_real tmp67;
	       fftw_real tmp98;
	       fftw_real tmp83;
	       fftw_real tmp102;
	       fftw_real tmp72;
	       fftw_real tmp99;
	       fftw_real tmp78;
	       fftw_real tmp101;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp64;
		    fftw_real tmp66;
		    fftw_real tmp63;
		    fftw_real tmp65;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp64 = X[4 * iostride];
		    tmp66 = Y[-5 * iostride];
		    tmp63 = c_re(W[3]);
		    tmp65 = c_im(W[3]);
		    tmp67 = (tmp63 * tmp64) - (tmp65 * tmp66);
		    tmp98 = (tmp65 * tmp64) + (tmp63 * tmp66);
	       }
	       {
		    fftw_real tmp80;
		    fftw_real tmp82;
		    fftw_real tmp79;
		    fftw_real tmp81;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp80 = X[iostride];
		    tmp82 = Y[-8 * iostride];
		    tmp79 = c_re(W[0]);
		    tmp81 = c_im(W[0]);
		    tmp83 = (tmp79 * tmp80) - (tmp81 * tmp82);
		    tmp102 = (tmp81 * tmp80) + (tmp79 * tmp82);
	       }
	       {
		    fftw_real tmp69;
		    fftw_real tmp71;
		    fftw_real tmp68;
		    fftw_real tmp70;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp69 = X[9 * iostride];
		    tmp71 = Y[0];
		    tmp68 = c_re(W[8]);
		    tmp70 = c_im(W[8]);
		    tmp72 = (tmp68 * tmp69) - (tmp70 * tmp71);
		    tmp99 = (tmp70 * tmp69) + (tmp68 * tmp71);
	       }
	       {
		    fftw_real tmp75;
		    fftw_real tmp77;
		    fftw_real tmp74;
		    fftw_real tmp76;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp75 = X[6 * iostride];
		    tmp77 = Y[-3 * iostride];
		    tmp74 = c_re(W[5]);
		    tmp76 = c_im(W[5]);
		    tmp78 = (tmp74 * tmp75) - (tmp76 * tmp77);
		    tmp101 = (tmp76 * tmp75) + (tmp74 * tmp77);
	       }
	       tmp73 = tmp67 - tmp72;
	       tmp84 = tmp78 - tmp83;
	       tmp85 = tmp73 + tmp84;
	       tmp91 = tmp67 + tmp72;
	       tmp92 = tmp78 + tmp83;
	       tmp93 = tmp91 + tmp92;
	       tmp100 = tmp98 + tmp99;
	       tmp103 = tmp101 + tmp102;
	       tmp128 = tmp100 + tmp103;
	       tmp121 = tmp98 - tmp99;
	       tmp122 = tmp101 - tmp102;
	       tmp142 = tmp121 + tmp122;
	  }
	  {
	       fftw_real tmp44;
	       fftw_real tmp105;
	       fftw_real tmp60;
	       fftw_real tmp109;
	       fftw_real tmp49;
	       fftw_real tmp106;
	       fftw_real tmp55;
	       fftw_real tmp108;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp41;
		    fftw_real tmp43;
		    fftw_real tmp40;
		    fftw_real tmp42;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp41 = X[2 * iostride];
		    tmp43 = Y[-7 * iostride];
		    tmp40 = c_re(W[1]);
		    tmp42 = c_im(W[1]);
		    tmp44 = (tmp40 * tmp41) - (tmp42 * tmp43);
		    tmp105 = (tmp42 * tmp41) + (tmp40 * tmp43);
	       }
	       {
		    fftw_real tmp57;
		    fftw_real tmp59;
		    fftw_real tmp56;
		    fftw_real tmp58;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp57 = X[3 * iostride];
		    tmp59 = Y[-6 * iostride];
		    tmp56 = c_re(W[2]);
		    tmp58 = c_im(W[2]);
		    tmp60 = (tmp56 * tmp57) - (tmp58 * tmp59);
		    tmp109 = (tmp58 * tmp57) + (tmp56 * tmp59);
	       }
	       {
		    fftw_real tmp46;
		    fftw_real tmp48;
		    fftw_real tmp45;
		    fftw_real tmp47;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp46 = X[7 * iostride];
		    tmp48 = Y[-2 * iostride];
		    tmp45 = c_re(W[6]);
		    tmp47 = c_im(W[6]);
		    tmp49 = (tmp45 * tmp46) - (tmp47 * tmp48);
		    tmp106 = (tmp47 * tmp46) + (tmp45 * tmp48);
	       }
	       {
		    fftw_real tmp52;
		    fftw_real tmp54;
		    fftw_real tmp51;
		    fftw_real tmp53;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp52 = X[8 * iostride];
		    tmp54 = Y[-iostride];
		    tmp51 = c_re(W[7]);
		    tmp53 = c_im(W[7]);
		    tmp55 = (tmp51 * tmp52) - (tmp53 * tmp54);
		    tmp108 = (tmp53 * tmp52) + (tmp51 * tmp54);
	       }
	       tmp50 = tmp44 - tmp49;
	       tmp61 = tmp55 - tmp60;
	       tmp62 = tmp50 + tmp61;
	       tmp88 = tmp44 + tmp49;
	       tmp89 = tmp55 + tmp60;
	       tmp90 = tmp88 + tmp89;
	       tmp107 = tmp105 + tmp106;
	       tmp110 = tmp108 + tmp109;
	       tmp127 = tmp107 + tmp110;
	       tmp118 = tmp105 - tmp106;
	       tmp119 = tmp108 - tmp109;
	       tmp141 = tmp118 + tmp119;
	  }
	  {
	       fftw_real tmp115;
	       fftw_real tmp86;
	       fftw_real tmp116;
	       fftw_real tmp124;
	       fftw_real tmp126;
	       fftw_real tmp120;
	       fftw_real tmp123;
	       fftw_real tmp125;
	       fftw_real tmp117;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp115 = K559016994 * (tmp62 - tmp85);
	       tmp86 = tmp62 + tmp85;
	       tmp116 = tmp39 - (K250000000 * tmp86);
	       tmp120 = tmp118 - tmp119;
	       tmp123 = tmp121 - tmp122;
	       tmp124 = (K951056516 * tmp120) + (K587785252 * tmp123);
	       tmp126 = (K951056516 * tmp123) - (K587785252 * tmp120);
	       Y[-5 * iostride] = tmp39 + tmp86;
	       tmp125 = tmp116 - tmp115;
	       Y[-7 * iostride] = tmp125 - tmp126;
	       X[3 * iostride] = tmp125 + tmp126;
	       tmp117 = tmp115 + tmp116;
	       Y[-9 * iostride] = tmp117 - tmp124;
	       X[iostride] = tmp117 + tmp124;
	  }
	  {
	       fftw_real tmp148;
	       fftw_real tmp143;
	       fftw_real tmp149;
	       fftw_real tmp147;
	       fftw_real tmp151;
	       fftw_real tmp145;
	       fftw_real tmp146;
	       fftw_real tmp152;
	       fftw_real tmp150;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp148 = K559016994 * (tmp141 - tmp142);
	       tmp143 = tmp141 + tmp142;
	       tmp149 = tmp144 - (K250000000 * tmp143);
	       tmp145 = tmp50 - tmp61;
	       tmp146 = tmp73 - tmp84;
	       tmp147 = (K951056516 * tmp145) + (K587785252 * tmp146);
	       tmp151 = (K587785252 * tmp145) - (K951056516 * tmp146);
	       X[5 * iostride] = -(tmp143 + tmp144);
	       tmp152 = tmp149 - tmp148;
	       X[7 * iostride] = tmp151 - tmp152;
	       Y[-3 * iostride] = tmp151 + tmp152;
	       tmp150 = tmp148 + tmp149;
	       X[9 * iostride] = -(tmp147 + tmp150);
	       Y[-iostride] = tmp150 - tmp147;
	  }
	  {
	       fftw_real tmp96;
	       fftw_real tmp94;
	       fftw_real tmp95;
	       fftw_real tmp112;
	       fftw_real tmp114;
	       fftw_real tmp104;
	       fftw_real tmp111;
	       fftw_real tmp113;
	       fftw_real tmp97;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp96 = K559016994 * (tmp90 - tmp93);
	       tmp94 = tmp90 + tmp93;
	       tmp95 = tmp87 - (K250000000 * tmp94);
	       tmp104 = tmp100 - tmp103;
	       tmp111 = tmp107 - tmp110;
	       tmp112 = (K951056516 * tmp104) - (K587785252 * tmp111);
	       tmp114 = (K951056516 * tmp111) + (K587785252 * tmp104);
	       X[0] = tmp87 + tmp94;
	       tmp113 = tmp96 + tmp95;
	       X[4 * iostride] = tmp113 - tmp114;
	       Y[-6 * iostride] = tmp113 + tmp114;
	       tmp97 = tmp95 - tmp96;
	       X[2 * iostride] = tmp97 - tmp112;
	       Y[-8 * iostride] = tmp97 + tmp112;
	  }
	  {
	       fftw_real tmp134;
	       fftw_real tmp129;
	       fftw_real tmp133;
	       fftw_real tmp138;
	       fftw_real tmp140;
	       fftw_real tmp136;
	       fftw_real tmp137;
	       fftw_real tmp139;
	       fftw_real tmp135;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp134 = K559016994 * (tmp127 - tmp128);
	       tmp129 = tmp127 + tmp128;
	       tmp133 = tmp132 - (K250000000 * tmp129);
	       tmp136 = tmp91 - tmp92;
	       tmp137 = tmp88 - tmp89;
	       tmp138 = (K951056516 * tmp136) - (K587785252 * tmp137);
	       tmp140 = (K951056516 * tmp137) + (K587785252 * tmp136);
	       Y[0] = tmp129 + tmp132;
	       tmp139 = tmp134 + tmp133;
	       X[6 * iostride] = -(tmp139 - tmp140);
	       Y[-4 * iostride] = tmp140 + tmp139;
	       tmp135 = tmp133 - tmp134;
	       X[8 * iostride] = -(tmp135 - tmp138);
	       Y[-2 * iostride] = tmp138 + tmp135;
	  }
     }
     if (i == m) {
	  fftw_real tmp1;
	  fftw_real tmp24;
	  fftw_real tmp8;
	  fftw_real tmp10;
	  fftw_real tmp25;
	  fftw_real tmp26;
	  fftw_real tmp14;
	  fftw_real tmp28;
	  fftw_real tmp23;
	  fftw_real tmp17;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = X[0];
	  tmp24 = X[5 * iostride];
	  {
	       fftw_real tmp2;
	       fftw_real tmp3;
	       fftw_real tmp4;
	       fftw_real tmp5;
	       fftw_real tmp6;
	       fftw_real tmp7;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp2 = X[4 * iostride];
	       tmp3 = X[6 * iostride];
	       tmp4 = tmp2 - tmp3;
	       tmp5 = X[8 * iostride];
	       tmp6 = X[2 * iostride];
	       tmp7 = tmp5 - tmp6;
	       tmp8 = tmp4 + tmp7;
	       tmp10 = K559016994 * (tmp4 - tmp7);
	       tmp25 = tmp2 + tmp3;
	       tmp26 = tmp5 + tmp6;
	  }
	  {
	       fftw_real tmp12;
	       fftw_real tmp13;
	       fftw_real tmp22;
	       fftw_real tmp15;
	       fftw_real tmp16;
	       fftw_real tmp21;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp12 = X[iostride];
	       tmp13 = X[9 * iostride];
	       tmp22 = tmp12 + tmp13;
	       tmp15 = X[3 * iostride];
	       tmp16 = X[7 * iostride];
	       tmp21 = tmp15 + tmp16;
	       tmp14 = tmp12 - tmp13;
	       tmp28 = K559016994 * (tmp22 + tmp21);
	       tmp23 = tmp21 - tmp22;
	       tmp17 = tmp15 - tmp16;
	  }
	  X[2 * iostride] = tmp1 + tmp8;
	  {
	       fftw_real tmp18;
	       fftw_real tmp20;
	       fftw_real tmp11;
	       fftw_real tmp19;
	       fftw_real tmp9;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp18 = (K587785252 * tmp14) - (K951056516 * tmp17);
	       tmp20 = (K951056516 * tmp14) + (K587785252 * tmp17);
	       tmp9 = tmp1 - (K250000000 * tmp8);
	       tmp11 = tmp9 - tmp10;
	       tmp19 = tmp10 + tmp9;
	       X[3 * iostride] = tmp11 - tmp18;
	       X[iostride] = tmp11 + tmp18;
	       X[4 * iostride] = tmp19 - tmp20;
	       X[0] = tmp19 + tmp20;
	  }
	  Y[-2 * iostride] = tmp23 - tmp24;
	  {
	       fftw_real tmp27;
	       fftw_real tmp32;
	       fftw_real tmp30;
	       fftw_real tmp31;
	       fftw_real tmp29;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp27 = (K951056516 * tmp25) + (K587785252 * tmp26);
	       tmp32 = (K951056516 * tmp26) - (K587785252 * tmp25);
	       tmp29 = (K250000000 * tmp23) + tmp24;
	       tmp30 = tmp28 + tmp29;
	       tmp31 = tmp29 - tmp28;
	       Y[0] = -(tmp27 + tmp30);
	       Y[-4 * iostride] = tmp27 - tmp30;
	       Y[-iostride] = tmp31 - tmp32;
	       Y[-3 * iostride] = tmp32 + tmp31;
	  }
     }
}

static const int twiddle_order[] =
{1, 2, 3, 4, 5, 6, 7, 8, 9};
fftw_codelet_desc fftw_hc2hc_forward_10_desc =
{
     "fftw_hc2hc_forward_10",
     (void (*)()) fftw_hc2hc_forward_10,
     10,
     FFTW_FORWARD,
     FFTW_HC2HC,
     223,
     9,
     twiddle_order,
};

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt