perlin.c

/*
perlin.c - Creates an image of Perlin (gradient) noise.
Default will write to a BMP. (Optionally writes it to a CSV file.)
Original Version: Alef Farah                     https://gist.github.com/a442/e82c5cbd29083558bd81
Cleaned up ver: Michaelangel007 aka MysticReddit https://gist.github.com/Michaelangel007/89bbb866cb048938df04
Enhancements:
  + Save perlin.bmp as default
  + Optional OpenMP support
  + Optional save perlin.CSV 
  + Cleaned up code for readability 
    + multi-column alignment
    + Synced x,y to i,j
  + Fixed array out-of-bounds in grad()
  + Fixed reversed arguments of lerp so it is consistent with "canonical" implementations
  + Fixed the backwards j,i and y,x so they are in normal order

Compilling:

  Without OpenMP
    gcc perlin.c -std=c99 -fopenmp -o perlin

  With OpenMP
    gcc perlin.c -std=c99 -fopenmp -o perlinmp -DOMP

OSX:
  Apple switched from gcc to llvm which doesn't support OpenMP as of May 2015.
  Install gcc 5.0
    curl -o http://iweb.dl.sourceforge.net/project/hpc/hpc/gcc/gcc-5.0-bin.tar.gz
    gunzip        gcc-5.0-bin.tar.gz
    sudo tar -xvf gcc-5.0-bin.tar    -C /.
    /usr/local/bin/gcc --version

Makefile: (You *must* preserve the tabs!)
- - - 8< - - -
all: perlin perlin_omp perlin_fix

clean:
	rm perlin perlin_omp perlin_fix

CC=/usr/local/bin/gcc
CFLAGS=-Wall -Wextra -std=c99
C_OMP=-DOMP=1
L_OMP=-fopenmp

perlin: perlin.c
	$(CC) $(CFLAGS)                $<          -o $@

perlin_omp: perlin.c
	$(CC) $(CFLAGS) $(C_OMP)       $< $(L_OMP) -o $@

perlin_fix: perlin.c
	$(CC) $(CFLAGS) -DFIXED_HASH=1 $< $(L_OMP) -o $@
- - - 8< - - -

Running:
  ./out [NUM_THREADS] [FNAME] [W] [H] [HASHSIZE] [GRIDSIZE] [SEED]

Details:
  Generate a WxH B&W image consisting of Perlin noise. Writes to a csv file
  called FNAME only if FNAME is passed and valid (e.g. of invalid in *nix:
  "/"). 8 different gradients are used, that is hardcoded. Assign one gradient
  for each point using a "hashtable" of size HASHSIZE, consisting of random
  numbers used to map points to gradients. Scale to GRIDSIZE before interpolating
  Use SEED in srand() (uint) and NUM_THREADS threads. See main() for defaults.

References:
  http://mrl.nyu.edu/~perlin/doc/oscar.html#noise
  http://mrl.nyu.edu/~perlin/noise/
  "Improving Noise"
    http://mrl.nyu.edu/~perlin/paper445.pdf
    Note: There is NO reference code nor reference data provided!
  http://webstaff.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf
  http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/perlin-noise-math-faq.html

Also of interest:
  http://devmag.org.za/2009/04/25/perlin-noise/
  http://stackoverflow.com/questions/8659351/2d-perlin-noise
*/

// Libs
  #include <stdlib.h>
  #include <stdio.h>
  #include <time.h>
  #include <stdint.h> // uint8_t, uint32_t
  #include <string.h> // strcat()

  #if OMP
    #include "omp.h"
  #else
    #define OMP 0     
  #endif

// Defines
  #ifndef   BMP
    #define BMP 1
  #endif

  #ifndef     USE_PERLIN_PERMUTATION
      #define USE_PERLIN_PERMUTATION 0 // use improved Perlin Noise permutations instead of random hashes
  #endif

// Types
  typedef short hash_t;

// Globals
  // Reference Perlin Noise Permutation Table (optional)
  static const hash_t PERMUTATION[] = { 
    151,160,137,91,90,15,
    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
    // Space-vs-Time tradeoff. Mirrored to prevent clamping & 0xFF
    151,160,137,91,90,15,
    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
  };

// Implementation

/*
   Original Perlin: 3x^2 - 2x^3
   Improved Perlin: 6x^5 - 15x^4 + 10x^3
*/
static inline double fade(double x) {
  return x*x*x*(x*(x*6 - 15) + 10);
}

/*
  Linear interpolation of a and b with a normalized percentage weight w
  Note: Also called 'mix'
  Note: You will want to compare speed and accuracy with your compiler,
        as a fused Multiply-Add might be available.
        https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation

        a       + w*(b-a)
        (1-w)*a + w*b
*/
static inline double lerp(double w,double a, double b) {
  return (1.0 - w)*a + w*b;
}

// The 8 unit gradients (center of a cube to its edges)
// Replaced gradient table with inlined add/sub operations
static inline double grad(const short hash, double x, double y)
{
  switch (hash & 7) {      // Gradient Table Size, sx*x + sy*y
    case 0: return -x - y; // {-1, -1},
    case 1: return  0 - y; // {-1,  0},
    case 2: return  x - y; // {-1,  1},
    case 3: return -x + 0; // { 0, -1},
    //      return  0 + 0; // { 0,  0} // not used
    case 4: return  x + 0; // { 0,  1},
    case 5: return -x + y; // { 1, -1},
    case 6: return  0 + y; // { 1,  0},
    case 7: return  x + y; // { 1,  1}
    default:return      0; // never happens
  }
}

/* 
  Improved Perlin 2D noise
  Bilinear interpolate on (x, y) using the influence of 4 gradients and
  a fade function as weight. 
  @param hashsize table's size; assumed to be a power of 2
  @param hashes "hashtable" used to pick gradients from the list
  @return interpolated noise value at <x, y> in range -1..+1 inclusive
*/
double noise2d(const size_t hashsize, const hash_t *hashtable, double x, double y) {
  // Get the coord of the top-left gradient of the grid (x,y) falls in
  const int    ix  = (int)x; 
  const int    iy  = (int)y;

  // Get the distance (x,y) is from it
  const double dx  = x - ix; // effective floor(): float -> int -> float
  const double dy  = y - iy;

  // Influence of gradient(ix,iy) (starting at the top-left one)
  // Do one necessary clamp to avoid expensive and unecessary % hashsize clamping
  const int    s   = ix & (hashsize-1); // Optimization: % hashsize
  const int    t   = iy & (hashsize-1);

  // Optimization: If we allow the table to be twice as big
  // then we can remove all remaining '% hashsize' clamping fluff
  const short  h0  = hashtable[(t  ) + hashtable[s  ]];
  const short  h1  = hashtable[(t  ) + hashtable[s+1]];
  const short  h2  = hashtable[(t+1) + hashtable[s  ]];
  const short  h3  = hashtable[(t+1) + hashtable[s+1]];

  const double g00 = grad( h0, dx  , dy   );
  const double g01 = grad( h1, dx-1, dy   );
  const double g10 = grad( h2, dx  , dy-1 );  
  const double g11 = grad( h3, dx-1, dy-1 );  

  const double fx  = fade(dx);
  const double fy  = fade(dy);

  /* Interpolate the influences using the blending function */
  const double y1  = lerp(fx, g00, g01); // top
  const double y2  = lerp(fx, g10, g11); // bot
  return             lerp(fy, y1 , y2 );
}

int main(int argc, char **argv) {
  char  filename[256] = "perlin";
  char *fname = filename;

  // filename
#if USE_PERLIN_PERMUTATION
    strcat( fname, "_perm" );
#endif
#if OMP
    strcat( fname, "_omp" );
#endif

  // extension
#if BMP
    strcat( fname, ".bmp" );
#else
    strcat( fname, ".csv" );
#endif

  size_t w = 800, h = 600, hashsize = 256, gridsize=32;
  int seed = !USE_PERLIN_PERMUTATION; // If using the Perlin permutations, don't perturb.

#if OMP
  unsigned int nthreads = omp_get_num_procs();
  if (argc >= 2 && atoi( argv[1]) > 0)
     nthreads =    atoi( argv[1]);
#endif

  if (argc >= 3) fname    = argv[2];
  if (argc >= 4) w        = atoi(argv[3]);
  if (argc >= 5) h        = atoi(argv[4]);
  if (argc >= 6) hashsize = atoi(argv[5]);
  if (argc >= 7) gridsize = atof(argv[6]);
  if (argc >= 8) seed     = atoi(argv[7]);

  if (!seed) {
       seed = (int)time(NULL);
  }
  srand((unsigned int)seed); // Note: srand(0) == srand(1)

  size_t i, j;

  // We allocate twice the space as an clamping optimization
  hash_t *hashes = (hash_t*)malloc(hashsize*2*sizeof(hash_t));

  // Generate the hashes
#if USE_PERLIN_PERMUTATION
  hashsize = 256;
  for (i = 0; i < hashsize; ++i) {
    hashes[i] = PERMUTATION[i];
  }
#else    
  for (i = 0; i < hashsize; ++i) {
    hashes[i] = i;
  }
#endif

  if( seed ) {
    /* Knuth shuffle (note that the hash values must be < hashsize for grad()) */
    for (int src = hashsize - 1; src > 0; --src) {
      int dst = rand() % src; // *ugh* Notoriously bad pseudo-random Linear Congruential Generator ...
      int tmp = hashes[dst];
                hashes[dst] = hashes[src];
                              hashes[src] = tmp;
    }
  }

  // Copy bottom 256 permutations to top 256 -- this is an optimization to avoid extensive clamping
  const int src =   0;
  const int dst = hashsize;
  for (i = 0; i < hashsize; i++ ) {
    hashes[ i + dst ] = hashes[ i + src ];
  }

  /* Allocate space for the img -- array of pointers */
  typedef uint32_t pix_t;
  pix_t **img = (pix_t**)malloc(h*sizeof(pix_t*));

#if OMP
  #pragma omp parallel for num_threads(nthreads) default(none) private(i) shared(h, w, img)
#endif
  for (j = 0; j < h; ++j) {
    img[j] = (pix_t*)malloc(w*sizeof(pix_t));
  }

  double x, y, n;
  const double r = 1.0 / (double)gridsize;

#if OMP
  #pragma omp parallel for collapse(2) num_threads(nthreads) default(none) private(i, j, x, y, n) shared(w, h, img, hashes, hashsize, gridsize)
#endif
  for (j = 0; j < h; ++j) { // *sigh* fixup reversed i,j and x,y
    for (i = 0; i < w; ++i) {
      x = i * r;
      y = j * r;

      n = noise2d(hashsize, hashes, x, y); // Optimization: remove constant GRADIENTS

      /* Scale the noise (which goes from -1 to 1) to [0, 255] */
      img[j][i] = (pix_t)((n + 1.0)*255.0)*0.5;
    }
  }

  free(hashes);

#if BMP
  printf( "Saving BMP: %s\n", fname );

  uint32_t headers[13];
  int      extrabytes = (w * 3) % 4;
  int      paddedsize = (w * 3 + extrabytes) * h;

  // https://msdn.microsoft.com/en-us/library/windows/desktop/dd183376%28v=vs.85%29.aspx
  // https://en.wikipedia.org/wiki/BMP_file_format
  headers[ 0] = paddedsize + 54; // bfSize (total)
  headers[ 1] = 0;               // bfReserved (both)
  headers[ 2] = 54;              // bfOffbits
  headers[ 3] = 40;              // biSize (sizeof BITMAPINFOHEADER)
  headers[ 4] = w;               // biWidth /sarcasm: Because someday BMPs will be larger then 64K x 64K, right!?
  headers[ 5] = h;               // biHeight
  headers[ 6] = (24 << 16) | 1;  // biPlanes=1 biBitCount=24
  headers[ 7] = 0;               // biCompression
  headers[ 8] = paddedsize;      // biSizeImage
  headers[ 9] = 0;               // biXPelsPerMeter
  headers[10] = 0;               // biYPelsPerMeter
  headers[11] = 0;               // biClrUsed
  headers[12] = 0;               // biClrImportant

  FILE *f = fopen( fname, "wb" );
  if( f ) {
    int n;
    fprintf(f, "BM"); // "MAGIC" 2 byte id as Bitmap

    for (n = 0; n <= 12; n++) {
       fprintf(f, "%c", (((uint32_t)headers[n]) >>  0) & 0xFF);
       fprintf(f, "%c", (((uint32_t)headers[n]) >>  8) & 0xFF);
       fprintf(f, "%c", (((uint32_t)headers[n]) >> 16) & 0xFF);
       fprintf(f, "%c", (((uint32_t)headers[n]) >> 24) & 0xFF);
    }

    // Stupid Windows BMP written upside down
    for( int y = h; --y >= 0; ) {
      const pix_t *scanline = &img[y][0];
      for( size_t x = 0; x < w; x++ ) {
        uint8_t pixel = *scanline++ & 0xFF;
        fprintf(f, "%c", pixel );
        fprintf(f, "%c", pixel );
        fprintf(f, "%c", pixel );
      }

      if (extrabytes) // See above - BMP lines must be of lengths divisible by 4
        for (n = 0; n < extrabytes; n++)
          fprintf(f, "%c", 0 );
    }
    fclose( f );
  }
#else
  printf( "Saving CSV: %s\n", fname );

  /* Print to file if requested and free */
  FILE *f = fopen(fname, "w");
  if( f ) {
    for (j = 0; j < h; ++j) { // *sigh* fixup reversed i,j and w,h
      for (i = 0; i < w; ++i) {
        fprintf(f, "%u,", img[j][i]);
      }
      fprintf(f, "\n");
    }
    fclose(f);
  }
#endif

  for (j = 0; j < h; j++)
    free( img[j] );
  free(img);

  return EXIT_SUCCESS;
}

I added some argument processing to my version, thought it might interest you if you're still maintaining yours.