Perlin Noise. Justin wrote this, I haven't actually looked at it yet... ANSI C

PB_perlin..h

#ifndef PB_PERLIN_H
#define PB_PERLIN_H

#include "PB_geom.h"
#include "PB_texture.h"
#include <stdlib.h>
#include <math.h>

typedef struct PerlinNoise
{
  Vec2 *grid; //size = gridWidth * gridHeight
  float octaveWidth;
  float octaveHeight;
  unsigned int seed;
  unsigned int gridWidth;
  unsigned int gridHeight;
  unsigned int numLevels;
}
PerlinNoise;

//float perlinNoise(const PerlinNoise *noise, const Vec2 sampleCoord, const Vec2 *offsets);

//PerlinNoise *myNoise = PerlinNoise_new(octaveWidth, octaveHeight, numOctaves, seed, gridWidth, gridHeight);

PerlinNoise *PerlinNoise_new
(
  const float octaveWidth, 
  const float octaveHeight,  
  const unsigned int seed, 
  const unsigned int gridWidth, 
  const unsigned int gridHeight,
  const unsigned int numLevels
);

void PerlinNoise_dispose(PerlinNoise *noise);

float perlinNoise(PerlinNoise *noise, const Vec2 sampleCoord, const Vec2 *offsets);

float sCurve(Vec2 r); // r is displacement

void perlinNoiseBatch
(
  const PerlinNoise *noise, 
  float *buffer, 
  const int width, 
  const int height, 
  const Vec2 *offsets
);

// noise function, must return same rand # if same # is given twice
// pretty much a random number function
float findNoise(float x, float y); 

// sets up the noise
float noise(float x, float y);

// interpolation function
float perlinInterpolate(Vec2 g, Vec2 r);

#endif

PB_perlin.c

#include "PB_perlin.h"

// function returning random number between [-1.0, 1.0]
// most common function found online for pseudo random noise generation
float findNoise(float x, float y)
{
  int number = (int)(x + y * 73); 
  int random_num;

  number ^= (number << 13); // bit shift and xor the number

  random_num = (number * (number * number * 60493 + 19990303) + 1376312589) 
    & 0x7fffffff; // and the number with this 8 bit

  return 1.0 - ((float) random_num / 1073741824.0);

}

float sCurve(Vec2 r)
{
  float x = 3 * (r.x * r.x) - 2 * (r.x * r.x * r.x); 
  float y = 3 * (r.y * r.y) - 2 * (r.y * r.y * r.y);

  return x * y;
}

float perlinInterpolate(Vec2 g, Vec2 displacement)
{
  float result = dot(g, displacement);

  result *= sCurve(displacement);

  return result;
}



PerlinNoise *PerlinNoise_new
  (
  const float octaveWidth, 
  const float octaveHeight,  
  const unsigned int seed, 
  const unsigned int gridWidth, 
  const unsigned int gridHeight,
  const unsigned int numLevels
  )

{
  int i = 0;
  int j = 0;
  PerlinNoise *perlinNoise = (PerlinNoise*)malloc(sizeof(PerlinNoise));

  perlinNoise->grid = (Vec2 *) malloc(sizeof(Vec2) * gridWidth * gridHeight);

  //fill in the noise for each index
  for ( i = 0; i < gridHeight; i++)
  {
    for (j = 0; j < gridWidth; j++)
    {
      perlinNoise->grid[j + (i * gridWidth)] = vec2(findNoise(seed + i, seed + j), findNoise(seed + i + 0.5, seed + j + 0.5));
      normalizeThis(perlinNoise->grid + (j + (i * gridWidth)));
    }
  }

  perlinNoise->octaveWidth = octaveWidth;
  perlinNoise->octaveHeight = octaveHeight;
  perlinNoise->seed = seed;
  perlinNoise->gridWidth = gridWidth;
  perlinNoise->gridHeight = gridHeight;
  perlinNoise->numLevels = numLevels;

  return perlinNoise;
}

void PerlinNoise_dispose(PerlinNoise *noise)
{
  free(noise->grid);
}

float perlinNoise(PerlinNoise *noise, const Vec2 sampleCoord, const Vec2 *offsets)
{
  int i = 0;
  int octaveFactor = 1 << (noise->numLevels - 1); // bit shift to get power of factor
  int power = 1;
  float result = 0;

  // adding noise levels, ie 4L0 + 2L1 + L2
  for (i = 0; i < noise->numLevels; ++i)
  {
    // wrapping ( point 100 % 101 equivalent is 1)
    float sampleX = fmodf(sampleCoord.x * power + offsets[i].x, noise->gridWidth * noise->octaveWidth);
    float sampleY = fmodf(sampleCoord.y * power + offsets[i].y, noise->gridHeight * noise->octaveHeight);

    float levelResult = 0;

    //grid[x + y * gridWidth] corresponds to (x * octaveWidth, y * octaveHeight)
    //the top left grid element relative to (sampleX, sampleY) is grid[((int) (sampleX / octaveWidth)) + ((int) (sampleY / octaveHeight)) * gridWidth]
    int indexTLX = (int) (sampleX / noise->octaveWidth);
    int indexTLY = (int) (sampleY / noise->octaveHeight);

    //displacement parameters
    float tx = 1.0f - fmodf(sampleX, noise->octaveWidth) / noise->octaveWidth;
    float ty = 1.0f - fmodf(sampleY, noise->octaveHeight) / noise->octaveHeight;

    // displacement vectors
    Vec2 disTL, disTR, disBL, disBR;

    // point vectors
    Vec2 TL, TR, BL, BR;

    TL = noise->grid[indexTLX + indexTLY * noise->gridWidth]; //top left
    TR = noise->grid[((indexTLX + 1) % noise->gridWidth) + indexTLY * noise->gridWidth]; // top right
    BL = noise->grid[indexTLX + ((indexTLY + 1) % noise->gridHeight) * noise->gridWidth];// bottom left
    BR = noise->grid[((indexTLX + 1) % noise->gridWidth) + ((indexTLY + 1) % noise->gridHeight) * noise->gridWidth];

    disTL = vec2(tx, ty);
    disTR = vec2(1.0f - tx, ty);
    disBL = vec2(tx, 1.0f - ty);
    disBR = vec2(1.0f - tx, 1.0f - ty);

    // get all point interpolations

    levelResult += perlinInterpolate(TL, disTL);
    levelResult += perlinInterpolate(TR, disTR);
    levelResult += perlinInterpolate(BL, disBL);
    levelResult += perlinInterpolate(BR, disBR);
    levelResult *= 0.5; // quick hack will fix later

    result += levelResult * octaveFactor;
    octaveFactor >>= 1; // divide by two
    power *= 2;
  }

  // divide by multiplying coeficients
  result /= (float) ((1 << noise->numLevels) - 1);

  return result;
}

// more optimized noise function
void perlinNoiseBatch(PerlinNoise * noise, float *buffer, const int width, const int height, const Vec2 *offsets)
{
  int i = 0;
  int octaveFactor = 1 << (noise->numLevels - 1);
  int power = 1;
  int x = 0;
  int y = 0;
  float divider = 1.0f / ((1 << noise->numLevels) - 1);

  for(x = 0; x < width; x++)
  {
    for(y = 0; y < height; y++)
    {
      buffer[x + y * width] = 0.0f;
    }
  }

  for (i = 0; i < noise->numLevels ;i++)
  {

    for (x = 0; x < width; x++)
    {
      float sampleX = fmod(x * power + offsets[i].x, noise->gridWidth * noise->octaveWidth);
      int floorX = (int)sampleX;
      int indexTLX = (int) (sampleX / noise->octaveWidth);
      float tx = 1.0f - fmod(sampleX, noise->octaveWidth) / noise->octaveWidth;

      for (y = 0; y < height; y++)
      {
        float result = 0.0f;
        float sampleY = fmod(y * power + offsets[i].y, noise->gridHeight * noise->octaveHeight);     
        int floorY = (int)sampleY;        
        int indexTLY = (int) (sampleY / noise->octaveHeight);        
        float ty = 1.0f - fmod(sampleY, noise->octaveWidth) / noise->octaveWidth;

        Vec2 disTL, disTR, disBL, disBR;

        Vec2 TL, TR, BL, BR;

        TL = noise->grid[indexTLX + indexTLY * noise->gridWidth];
        TR = noise->grid[((indexTLX + 1) % noise->gridWidth) + indexTLY * noise->gridWidth];
        BL = noise->grid[indexTLX + ((indexTLY + 1) % noise->gridHeight) * noise->gridWidth];
        BR = noise->grid[((indexTLX + 1) % noise->gridWidth) + ((indexTLY + 1) % noise->gridHeight) * noise->gridWidth];

        disTL = vec2(tx, ty);
        disTR = vec2(1.0f - tx, ty);
        disBL = vec2(tx, 1.0f - ty);
        disBR = vec2(1.0f - tx, 1.0f - ty);

        result += perlinInterpolate(TL, disTL);
        result += perlinInterpolate(TR, disTR);
        result += perlinInterpolate(BL, disBL);
        result += perlinInterpolate(BR, disBR);
        result *= 0.5; // quick hack that will be fixed later

        buffer[x + y * width] += result * octaveFactor;

      }
    }

    octaveFactor >>= 1; 
    power *= 2;
  }

  for(x = 0; x < width; x++)
  {
    for(y = 0; y < height; y++)
    {
      buffer[x + y * width] *= divider;
    }
  }

}