A simple Perlin Noise generator using OpenMP for parallelization

perlin.c

/*
 * perlin.c - create an image of Perlin (gradient) noise
 *
 * Compilling:
 * gcc perlin.c -std=c99 -O3 -Wall -Wextra -Werror -fopenmp -lgomp
 *
 * Doc:
 * ./a.out --help
 *
 * Thanks to:
 *  https://github.com/Michaelangel007 for bugfix.
 */

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <argp.h>
#include <omp.h>

/*
 * begin argp (argument processing) stuff
 */

#define ARGP_FLAGS 0
#define ARGP_INDEX 0
#define ARGP_MAX_ARGS 0
static char const ARGP_DOC[] = "Calculates some Perlin (gradient) noise, "
  "optionally outputting a .pgm to stdout.";
static char const ARGP_DOCA[] = "[NUM_THREADS] [OUTPUT] [WIDTH] [HEIGHT] "
  "[HSIZE] [GSIZE] [SEED]";
static struct argp_option const ARGP_OPT[] = {
  {"num_threads", 't', "NUM_THREADS", 0, "Number of threads to be used - "
    "default is the number of logical cores.", 0},
  {"output", 'o', NULL, OPTION_ARG_OPTIONAL, "Output a .pgm to stdout", 0},
  {"width", 'w', "WIDTH", 0, "Image width in pixels. Default is 800.", 0},
  {"height", 'h', "HEIGHT", 0, "Image height in pixels. Default is 600.", 0},
  {"hsize", 'z', "HSIZE", 0, "Size of the \"hash table\" of random numbers "
    "used to map gradients to points (see Ken Perlin's original algorithm). "
    "Default is 256 elements.", 0},
  {"gsize", 'g', "GSIZE", 0, "Scale to GSIZE before interpolating. Default "
    "is 32.", 0},
  {"seed", 's', "SEED", 0, "Seed srand with SEED. Default is time(NULL).", 0},
  { 0 }
};

struct argp_arguments {
  char *input[ARGP_MAX_ARGS];
  int output;
  unsigned num_threads, seed;
  size_t width, height, hsize;
  double gsize;
};

static error_t
argp_parse_options(int key, char *arg, struct argp_state *state)
{
  /* FIXME assumes values > 0 and sizeof(size_t) <= sizeof(long) */
  struct argp_arguments *arguments = (struct argp_arguments *)(state->input);
  switch(key) {
    case 'o': arguments->output = 1; break;
    case 't': arguments->num_threads = atol(arg); break;
    case 'w': arguments->width = atol(arg); break;
    case 'h': arguments->height = atol(arg); break;
    case 'z': arguments->hsize = atol(arg); break;
    case 'g': arguments->gsize = atof(arg); break;
    case 's': arguments->seed = atol(arg); break;
    case ARGP_KEY_ARG: argp_usage(state); break;
    case ARGP_KEY_END: return 0;
    default: return ARGP_ERR_UNKNOWN;
  }
  return 0;
}

/*
 * end argp stuff
 */

/* The fade function, 6x^5 - 15x^4 + 10x^3 */
static inline double
fade(double x)
{
  return x*x*x*(x*(x*6 - 15) + 10);
}

/* dot product between the gradient [gy, gx] and [y, x] */
static inline double
dot(int *g, double y, double x)
{
  return g[0]*y + g[1]*x;
}

/* Linear interpolation of a and b with weight w */
static inline double
lerp(double a, double b, double w)
{
  return (1 - w)*a + w*b;
}

/* Given the "hashtable" of size hsize, determine the gradient for (i, j) */
static inline int *
grad(size_t *hashes, size_t hsize, int grads[][2], int i, int j)
{
  return grads[hashes[(i % hsize + hashes[j % hsize]) % hsize] % 8];
}

/*
 * Bilinear interpol on (y, x) using the influence of 4 gradients and a fade
 * function as weight. Grads is a list of gradients, hashes the "hashtable"
 * used to pick gradients from the list and hsize the table's size. Returns the
 * interpolated noise value at (y, x) in [-1, 1].
 */
static double
interpol(size_t *hashes, size_t hsize, int grads[][2], double y, double x)
{
  /* Get the coord of the top-left gradient of the grid (y, x) falls in */
  int j = (int)x;
  int i = (int)y;
  /* Get the distance (y, x) is from it */
  double dx = x - j;
  double dy = y - i;
  /* Influence of (g)radient(i)(j) (starting at the top-left one) */
  double g00 = dot(grad(hashes, hsize, grads, i, j), dy, dx);
  double g01 = dot(grad(hashes, hsize, grads, i, j + 1), dy, dx - 1);
  double g10 = dot(grad(hashes, hsize, grads, i + 1, j), dy - 1, dx);
  double g11 = dot(grad(hashes, hsize, grads, i + 1, j + 1), dy - 1, dx - 1);
  /* Interpolate the influences using the blending function */
  /* Linear interpol the top 2 */
  double lt = lerp(g00, g01, fade(dx));
  /* Linear interpol the bottom 2 */
  double lb = lerp(g10, g11, fade(dx));
  /* Linear interpol lb lt, completing the bilienear interpol */
  return lerp(lt, lb, fade(dy));
}

static inline void
output_pgm(unsigned const *img, size_t w, size_t h, unsigned max_color)
{
  printf("P2 %zu %zu %u ", w, h, max_color);
  for (size_t i = 0; i < w * h; ++i)
    printf("%u ", img[i]);
}

static inline void
swap(size_t *a, size_t *b)
{
  size_t tmp = *a;
  *a = * b;
  *b = tmp;
}

int
main(int argc, char **argv)
{
  /* Default and user-defined parameters */
  struct argp_arguments args;
  args.output = 0;
  args.width = 800;
  args.height = 600;
  args.hsize = 256;
  args.gsize = 32;
  args.seed = time(NULL);
  args.num_threads = omp_get_num_procs();
  struct argp argp = {
    ARGP_OPT, argp_parse_options, ARGP_DOCA, ARGP_DOC, 0, 0, 0
  };
  if (argp_parse(&argp, argc, argv, ARGP_FLAGS, ARGP_INDEX, &args) ==
      ARGP_KEY_ERROR) {
    fprintf(stderr, "%s, error while parsing parameters\n", argv[0]);
    return EXIT_FAILURE;
  }
  srand(args.seed);
  /* Create 8 unit gradients (center of a cube to its edges) */
  int grads[8][2] = {
    {-1, -1}, {-1, 0}, {-1, 1}, {0, -1}, {0, 1}, {1, -1}, {1, 0}, {1, 1}
  };
  /*
   * Create a "hash table" with random numbers that a function will use to
   * map coordinates to a gradient pseudo-randomly.
   */
  size_t *hashes = malloc(args.hsize * sizeof(*hashes));
  for (size_t i = 0; i < args.hsize; ++i)
    hashes[i] = i;
  /* Knuth shuffle (note that the hash values must be < hsize for grad()) */
  for (size_t i = args.hsize - 1; i > 0; --i)
    swap(hashes + i, hashes + (size_t)(rand() % i));
  /* Allocate space for img */
  double inv_gsize = 1 / args.gsize;
  unsigned *img = malloc(args.width * args.height * sizeof(*img));
  /*
   * With OpenMP we don't want to collapse the loops because of the row-level
   * calculations; with a GPU we might have to think of a different strategy
   */
#pragma omp parallel for num_threads(args.num_threads) default(none)\
shared(img, args, hashes, grads, inv_gsize) schedule(static)
  for (size_t i = 0; i < args.height; ++i) {
    size_t row = i * args.width;
    double y = i * inv_gsize;
    for (size_t j = 0; j < args.width; ++j)
      /* Scale the noise from [-1, 1] to [0, 255] */
      img[row + j] = ((interpol(hashes, args.hsize, grads, y, j * inv_gsize) +
            1) * 255) * 0.5;
  }
  free(hashes);
  if (args.output)
    output_pgm(img, args.width, args.height, 255);
  free(img);
  return EXIT_SUCCESS;
}