Sorry your browser is not supported!

` Green Gandalf's Perlin style noise function - Version 3
` Created 25 June 2007, modified 17 February 2008.

`   Uses suggestions from following website
`      http://www.mandelbrot-dazibao.com/Perlin/Perlin1.htm

`   basic idea works with simple bilinear interpolation
`    - but has visible interior seams
`   doesn't work well at all with Perlin style interpolation

sync on: sync rate 0: sync
set display mode 800, 600, 32

randomize 140549 ` arbitrary fixed number for reproducibility

autocam off
position camera 0, 50, -300
point camera 0, 0, 0

create bitmap 1, 512, 512

global nOct = 8
global nGrid
global twoPiInv as float
twoPiInv = 0.159154943

nGrid = 2^nOct

dim weight#(nOct)
dim interpWeight(3) as float
dim tangent(3) as float

dim rawNoise(511, 511) as float

p# = 2.0
for i=0 to nOct
  weight#(i) = 1.0/p#^i
next i

dim a(nGrid, nGrid, nOct) as float
dim b(nGrid, nGrid, nOct) as float

for oct = 0 to nOct
  for i = 0 to nGrid
    for j = 0 to nGrid
      a(i, j, oct) = rnd(100) * 0.01  ` random value in range 0 to 1
      b(i, j, oct) = rnd(100) * 0.01  ` random value in range 0 to 1
                                             ` not seamless yet
    next j
  next i
next oct

` calculate raw noise and find scale factors for
` rescaling noise values to byte range
minNoise# = 1000  ` arbitrary large number out of range
maxNoise# = -1000 ` arbitrary small number out of range
for x = 0 to 511
  for y = 0 to 511
    u# = x * 0.001953125 ` value in range 0 to 1
    v# = y * 0.001953125
    rawNoise(x, y) = noise(u#, v#)
    if minNoise# > rawNoise(x, y)
      minNoise# = rawNoise(x, y)
    else
      if maxNoise# < rawNoise(x, y) then maxNoise# = rawNoise(x, y)
    endif
  next y
next x

f1# = 255.0/(maxNoise# - minNoise#)
f2# = minNoise# * 255.0/(maxNoise# - minNoise#)

lock pixels
  for x = 0 to 511
    for y = 0 to 511
      ` convert raw noise to byte range 0 - 255
      c = rawNoise(x, y) * f1# - f2#
      ` just in case
      if c<0 then c=0
      if c>255 then c=255
      dot x, y, rgb(0, c, 0) ` a nice green colour :)
    next y
  next x
unlock pixels

copy bitmap 1, 0, 0, 511, 511,  0, 144, 44, 655, 555

set current bitmap 0

repeat
  text 20, 20, "All done!"
  sync
until spacekey()

set current bitmap 1
get image 1, 0, 0, 512, 512
save image "test v3.png", 1
end

function noiseBase(oct, x as float, y as float)
  ` x and y both in range 0 to 1
  ` find correct "tile" to use
  tGrid = 2^oct
  x# = x * tGrid
  y# = y * tGrid
  i = floor(x#)
  if i >= tGrid then i = tGrid - 1 ` just in case
  j = floor(y#)
  if j >= tGrid then j = tGrid - 1 ` just in case
  x# = x# - i  ` number in range 0 to 1
  y# = y# - j
  ` calculate interpolation weights
  ` Perlin weights - image not right at all WHY?
  remstart
  ix1# = x# * x# * (3.0 - 2.0 * x#)
  ix0# = 1.0 - ix1#
  iy1# = y# * y# * (3.0 - 2.0 * y#)
  iy0# = 1.0 - iy1#
  remend
  ` simple bilinear weights - should show slight seams
  ` because derivatives not matched correctly
`  remstart
  ix0# = 1.0 - x#
  ix1# = x#
  iy0# = 1.0 - y#
  iy1# = y#
`  remend
  ` simple sine weights (a variation of the Perlin idea)
  remstart
  ix1# = x# - sin(360.0 * x#) * twoPiInv
  ix0# = 1.0 - ix1#
  iy1# = y# - sin(360.0 * y#) * twoPiInv
  iy0# = 1.0 - iy1#
  remend
  interpWeight(0) = ix0# * iy0#
  interpWeight(1) = ix0# * iy1#
  interpWeight(2) = ix1# * iy0#
  interpWeight(3) = ix1# * iy1#
  ` calculate tangents at four corners of tile
  tangent(0) = a(i, j, oct) * x# + b(i, j, oct) * y#
  tangent(1) = a(i, j + 1, oct) * x# + b(i, j + 1, oct) * (y# - 1.0)
  tangent(2) = a(i + 1 , j, oct) * (x# - 1.0) + b(i + 1, j, oct) * y#
  tangent(3) = a(i + 1, j + 1, oct) * (x# - 1.0) + b(i + 1, j + 1, oct) * (y# - 1.0)
  ` calculate interpolated noise value
  result# = 0.0
  for corner = 0 to 3
    result# = result# + interpWeight(corner) * tangent(corner)
  next corner
endfunction result#

function noise(x as float, y as float)
  result# = 0.0
`  result# = noiseBase(1, x, y)
  for oct = 0 to nOct
    result# = noiseBase(oct, x, y) * weight#(oct) + result#
  next oct
endfunction result#

  ` calculate interpolation weights
  ` Perlin weights - image not right at all WHY?
  remstart
  ix1# = x# * x# * (3.0 - 2.0 * x#)
  ix0# = 1.0 - ix1#
  iy1# = y# * y# * (3.0 - 2.0 * y#)
  iy0# = 1.0 - iy1#
  remend
  ` simple bilinear weights - should show slight seams
  ` because derivatives not matched correctly
`  remstart
  ix0# = 1.0 - x#
  ix1# = x#
  iy0# = 1.0 - y#
  iy1# = y#
`  remend

Load DLL "user32.dll"   , 1
displayX = Call DLL( 1 , "GetSystemMetrics" , 0 )
displayY = Call DLL( 1 , "GetSystemMetrics" , 1 )
Delete DLL 1

Set Display Mode displayX, displayY, 32, 1
Sync On

// Seed
Randomize 1337

// You can manually specify the resolution for higher detail noise
`startingRes     = 512
`totalOctaves    = 9

totalOctaves    = 8
startingRes     = 2 ^ totalOctaves

benchmark1Start = Timer()

// Generate noise for all octaves
for octaveID = 1 to totalOctaves
    memblockID  = octaveID
    octaveRes   = startingRes / 2 ^ (octaveID - 1)
    Make Memblock memblockID, octaveRes * octaveRes * 4 // A single float
    for x = 0 to octaveRes - 1
        for y = 0 to octaveRes - 1
            Write Memblock Float memblockID, ( x + y * octaveRes ) * 4, ( Rnd(2000) - 1000 ) * 0.001
        next
    next
next

benchmark1elapsed   = Timer() - benchmark1Start

benchmark2Start = Timer()

// add all layers together
memblockID  = totalOctaves + 1

Make Memblock memblockID, 12 + startingRes * startingRes * 4

Write Memblock DWord memblockID, 0, startingRes
Write Memblock DWord memblockID, 4, startingRes
Write Memblock DWord memblockID, 8, 32

minHeight#      = -1.9
maxHeight#      =  1.9
minInv#         = -minHeight#
heightRange#    =  maxHeight# - minHeight#
_255DivRange#   =  255.0 / heightRange#
sResReciprocal# =    1.0 / (startingRes)

for x = 0 to startingRes - 1
    for y = 0 to startingRes - 1
        x#          = x * sResReciprocal#
        y#          = y * sResReciprocal#
        height#     = 0.0
        Strength#   = 1.0
        `height# = Memblock Float( 1, ( x + y * startingRes ) * 4 )

        for OctaveID = totalOctaves to 1 step -1
        `for OctaveID = 7 to 7 step -1
            octaveRes = startingRes / 2 ^ (octaveID - 1)
            // Store the float pixel we are over on the current octave
            memblockX# = x# * octaveRes
            memblockY# = y# * octaveRes
            // Store the int version for pixel sampling
            memblockX  = memblockX#
            memblockY  = memblockY#
            // Get the local offset
            memblockX# = memblockX# mod 1.0
            memblockY# = memblockY# mod 1.0

            // BICUBIC ///////////
            `REMSTART
            sample1#  = Sample( OctaveID, memblockX - 1, memblockY - 1, octaveRes )
            sample2#  = Sample( OctaveID, memblockX    , memblockY - 1, octaveRes )
            sample3#  = Sample( OctaveID, memblockX + 1, memblockY - 1, octaveRes )
            sample4#  = Sample( OctaveID, memblockX + 2, memblockY - 1, octaveRes )
            mHeight1# = CubicInterpolate( sample1#, sample2#, sample3#, sample4#, memblockX# )

            sample1#  = Sample( OctaveID, memblockX - 1, memblockY    , octaveRes )
            sample2#  = Sample( OctaveID, memblockX    , memblockY    , octaveRes )
            sample3#  = Sample( OctaveID, memblockX + 1, memblockY    , octaveRes )
            sample4#  = Sample( OctaveID, memblockX + 2, memblockY    , octaveRes )
            mHeight2# = CubicInterpolate( sample1#, sample2#, sample3#, sample4#, memblockX# )

            sample1#  = Sample( OctaveID, memblockX - 1, memblockY + 1, octaveRes )
            sample2#  = Sample( OctaveID, memblockX    , memblockY + 1, octaveRes )
            sample3#  = Sample( OctaveID, memblockX + 1, memblockY + 1, octaveRes )
            sample4#  = Sample( OctaveID, memblockX + 2, memblockY + 1, octaveRes )
            mHeight3# = CubicInterpolate( sample1#, sample2#, sample3#, sample4#, memblockX# )

            sample1#  = Sample( OctaveID, memblockX - 1, memblockY + 2, octaveRes )
            sample2#  = Sample( OctaveID, memblockX    , memblockY + 2, octaveRes )
            sample3#  = Sample( OctaveID, memblockX + 1, memblockY + 2, octaveRes )
            sample4#  = Sample( OctaveID, memblockX + 2, memblockY + 2, octaveRes )
            mHeight4# = CubicInterpolate( sample1#, sample2#, sample3#, sample4#, memblockX# )

            mHeight#  = CubicInterpolate( mHeight1#, mHeight2#, mHeight3#, mHeight4#, memblockY# )
            `REMEND
            //////////////////////

            // BILINEAR //////////
            REMSTART
            sample1#  = Sample( OctaveID, memblockX    , memblockY, octaveRes )
            sample2#  = Sample( OctaveID, memblockX + 1, memblockY, octaveRes )
            mHeightX# = Linear_Interpolate( sample1#, sample2#, memblockX# )

            sample1#  = Sample( OctaveID, memblockX    , memblockY + 1, octaveRes)
            sample2#  = Sample( OctaveID, memblockX + 1, memblockY + 1, octaveRes)
            mHeightY# = Linear_Interpolate( sample1#, sample2#, memblockX# )

            mHeight#  = Linear_Interpolate( mHeightX#, mHeightY#, memblockY# )
            REMEND
            //////////////////////

            height#   = height# + mHeight# * Strength#
            // Make the next octave affect the output by half as much as this one
            Strength# = Strength# * 0.5
        next

        height  = (height# + minInv#) * _255DivRange#
        if height > 255 then height = 255
        if height <   0 then height = 0
        Write Memblock Byte memblockID, 12 + ( x + y * startingRes ) * 4    , height // B
        Write Memblock Byte memblockID, 12 + ( x + y * startingRes ) * 4 + 1, height // G
        Write Memblock Byte memblockID, 12 + ( x + y * startingRes ) * 4 + 2, height // R
        Write Memblock Byte memblockID, 12 + ( x + y * startingRes ) * 4 + 3, 255    // A
    next
next

benchmark2elapsed   = Timer() - benchmark2Start

// Create image
Make Image From Memblock 1, memblockID

// Loop to output result
do
    CLS

    Paste Image 1, 64, 64

    Text 5,  5, "Time to generate noise: " + Str$(benchmark1elapsed) + "ms"
    Text 5, 15, "Time to add octaves: "    + Str$(benchmark2elapsed) + "ms"

    Sync
loop

END

Function Sample( memblockID, x, y, memblockRes )

    // Wrap
    if x < 0 then x = x + memblockRes
    if y < 0 then y = y + memblockRes

    if x > memblockRes-1 then x = x - memblockRes
    if y > memblockRes-1 then y = y - memblockRes

    returnValue# = Memblock Float( memblockID, ( x + y * memblockRes ) * 4 )

Endfunction returnValue#

Function Linear_Interpolate( a as float, b as float, x as float )

    returnValue# =  a * (1.0 - x) + b * x

Endfunction returnValue#


Function CubicInterpolate( aMinOne as float, a as float, b as float, bAddOne as float, acrossAB as float )

    P as float
    Q as float
    R as float
    S as float

    P = (bAddOne - b) - (aMinOne - a)
    Q = (aMinOne - a) - P
    R = b - aMinOne
    S = a

    returnValue# =  P * acrossAB ^ 3 + Q * acrossAB ^ 2 + R * acrossAB + S

EndFunction returnValue#

#include "DarkGDK.h"
#include "Windows.h"
#include "Winuser.h"
#include <string>
#include "Math.h"

int displayX = GetSystemMetrics(SM_CXSCREEN);
int displayY = GetSystemMetrics(SM_CYSCREEN);

#define startingRes  512
#define totalOctaves 8

float* noiseMap[totalOctaves+1];

char globalString[256];

//#define INTERPOLATION_LINEAR
//#define INTERPOLATION_COSINE
#define INTERPOLATION_CUBIC

// Prototypes
float Sample( int octaveID, int x, int y, int octaveRes );
float LinearInterpolate( float a, float b, float x );
float CubicInterpolate( float aMinOne, float a, float b, float bAddOne, float acrossAB );
float CosineInterpolate( float a, float b, float x );

void DarkGDK ( void )
{
	// Set display
	dbSetWindowPosition( 0, 0 );
	dbSetDisplayMode( displayX, displayY, 32 );
	dbSetWindowLayout( 0, 0, 0 );
	
	// Init
	dbSyncOn();
	
	dbRandomize(1337);
	
	int benchmark1Start = dbTimer();
	
	// Generate noise for all octaves
	for( int octaveID = 1; octaveID <= totalOctaves; octaveID ++ )
	{
		int octaveRes   = startingRes / pow( 2.0f, (float) (octaveID - 1) );
		noiseMap[octaveID] = new float[octaveRes * octaveRes];
		int resLoopMax = octaveRes - 1;
		for( int x = 0; x <= resLoopMax; x ++ )
		{
			for( int y = 0; y <= resLoopMax; y ++ )
			{
				noiseMap[octaveID][x+y*octaveRes] = ( ((float) dbRnd(2000)) - 1000.0f ) * 0.001;
			}
		}
	}
	
	int benchmark1Elapsed = dbTimer() - benchmark1Start;
	
	int benchmark2Start = dbTimer();
	
	int memblockID  = 1;
	
	dbMakeMemblock( memblockID, 12 + startingRes * startingRes * 4 );
	
	dbWriteMemblockDword( memblockID, 0, startingRes );
	dbWriteMemblockDword( memblockID, 4, startingRes );
	dbWriteMemblockDword( memblockID, 8, 32			 );
	
	float minHeight			= -1.9;
	float maxHeight			=  1.9;
	float minInv			= -minHeight;
	float heightRange		=  maxHeight - minHeight;
	float _255DivRange		=  255.0 / heightRange;
	float sResReciprocal	=    1.0 / (startingRes);
	
	int resLoopMax = startingRes - 1;
	
	for( int x = 0; x <= resLoopMax; x ++ )
	{
		for( int y = 0; y <= resLoopMax; y ++ )
		{
			float fX		= ((float) x) * sResReciprocal;
			float fY		= ((float) y) * sResReciprocal;
			float height	= 0.0;
			float strength	= 1.0;
			float mHeight;
			
			for ( int octaveID = totalOctaves; octaveID >= 1; octaveID -- )
			{
				int   octaveRes = startingRes / pow( 2.0f, (float) (octaveID - 1) );
				float octavefX = fX * (float) octaveRes;
				float octavefY = fY * (float) octaveRes;
				int	  octaveX  = (int) octavefX;
				int	  octaveY  = (int) octavefY;
				octavefX = fmod( octavefX, 1.0f);
				octavefY = fmod( octavefY, 1.0f);
			
				// BICUBIC ///////////
				#ifdef INTERPOLATION_CUBIC
					float sample1  = Sample( octaveID, octaveX - 1, octaveY - 1, octaveRes );
					float sample2  = Sample( octaveID, octaveX    , octaveY - 1, octaveRes );
					float sample3  = Sample( octaveID, octaveX + 1, octaveY - 1, octaveRes );
					float sample4  = Sample( octaveID, octaveX + 2, octaveY - 1, octaveRes );
					float mHeight1 = CubicInterpolate( sample1, sample2, sample3, sample4, octavefX );

						  sample1  = Sample( octaveID, octaveX - 1, octaveY    , octaveRes );
						  sample2  = Sample( octaveID, octaveX    , octaveY    , octaveRes );
						  sample3  = Sample( octaveID, octaveX + 1, octaveY    , octaveRes );
						  sample4  = Sample( octaveID, octaveX + 2, octaveY    , octaveRes );
					float mHeight2 = CubicInterpolate( sample1, sample2, sample3, sample4, octavefX );

						  sample1  = Sample( octaveID, octaveX - 1, octaveY + 1, octaveRes );
						  sample2  = Sample( octaveID, octaveX    , octaveY + 1, octaveRes );
						  sample3  = Sample( octaveID, octaveX + 1, octaveY + 1, octaveRes );
						  sample4  = Sample( octaveID, octaveX + 2, octaveY + 1, octaveRes );
					float mHeight3 = CubicInterpolate( sample1, sample2, sample3, sample4, octavefX );

						  sample1  = Sample( octaveID, octaveX - 1, octaveY + 2, octaveRes );
						  sample2  = Sample( octaveID, octaveX    , octaveY + 2, octaveRes );
						  sample3  = Sample( octaveID, octaveX + 1, octaveY + 2, octaveRes );
						  sample4  = Sample( octaveID, octaveX + 2, octaveY + 2, octaveRes );
					float mHeight4 = CubicInterpolate( sample1, sample2, sample3, sample4, octavefX );

					mHeight  = CubicInterpolate( mHeight1, mHeight2, mHeight3, mHeight4, octavefY );
				#endif
				//////////////////////
				// COSINE ////////////
				#ifdef INTERPOLATION_COSINE
					float sample1  = Sample( octaveID, octaveX    , octaveY, octaveRes );
					float sample2  = Sample( octaveID, octaveX + 1, octaveY, octaveRes );
					float mHeightX = CosineInterpolate( sample1, sample2, octavefX );
					
						  sample1  = Sample( octaveID, octaveX    , octaveY + 1, octaveRes);
						  sample2  = Sample( octaveID, octaveX + 1, octaveY + 1, octaveRes);
					float mHeightY = CosineInterpolate( sample1, sample2, octavefX );
					
					mHeight = CosineInterpolate( mHeightX, mHeightY, octavefY );
				#endif
				/////////////////////
				// BILINEAR /////////
				#ifdef INTERPOLATION_LINEAR
					float sample1  = Sample( octaveID, octaveX    , octaveY, octaveRes );
					float sample2  = Sample( octaveID, octaveX + 1, octaveY, octaveRes );
					float mHeightX = LinearInterpolate( sample1, sample2, octavefX );
					
						  sample1  = Sample( octaveID, octaveX    , octaveY + 1, octaveRes);
						  sample2  = Sample( octaveID, octaveX + 1, octaveY + 1, octaveRes);
					float mHeightY = LinearInterpolate( sample1, sample2, octavefX );
					
					mHeight = LinearInterpolate( mHeightX, mHeightY, octavefY );
				#endif
				/////////////////////
				
				height   = height + mHeight * strength;
				strength = strength * 0.5f;
			}
			height   = (height + minInv) * _255DivRange;
			if (height > 255)
				height = 255;
			if (height <   0)
				height = 0;
			int iHeight = (int) height;
			dbWriteMemblockByte( memblockID, 12 + ( x + y * startingRes ) * 4    , iHeight ); // B
			dbWriteMemblockByte( memblockID, 12 + ( x + y * startingRes ) * 4 + 1, iHeight ); // G
			dbWriteMemblockByte( memblockID, 12 + ( x + y * startingRes ) * 4 + 2, iHeight ); // R
			dbWriteMemblockByte( memblockID, 12 + ( x + y * startingRes ) * 4 + 3, 255    ); // A
		}
	}
	
	dbMakeImageFromMemblock( 1, memblockID );
	
	int benchmark2Elapsed = dbTimer() - benchmark1Start;
	
	// Output image and info
	while ( LoopGDK() )
	{
		sprintf( globalString, "Time to generate noise: %ims", benchmark1Elapsed );
		dbText( 5, 5, globalString );
		sprintf( globalString, "Time to add octaves: %ims", benchmark2Elapsed );
		dbText( 5, 15, globalString );
		
		dbPasteImage( 1, 64, 64 );
		
		dbSync();
	}
	
	return;
}

float Sample( int octaveID, int x, int y, int octaveRes )
{
	// Wrap
	if (x < 0)
		x += octaveRes;
	if (y < 0)
		y += octaveRes;
	
	if (x > octaveRes-1)
		x -= octaveRes;
	if (y > octaveRes-1)
		y -= octaveRes;
	
	return noiseMap[octaveID][x+y*octaveRes];
}

float LinearInterpolate( float a, float b, float x )
{
	return a * (1.0 - x) + b * x;
}
float CosineInterpolate( float a, float b, float x )
{
    float f = (1.0 - cos(x * 3.141592654)) * 0.5;
    
	return a * (1.0 - f) + b * f;
}
float CubicInterpolate( float aMinOne, float a, float b, float bAddOne, float acrossAB )
{
	float P, Q, R, S;
	
	P = (bAddOne - b) - (aMinOne - a);
	Q = (aMinOne - a) - P;
	R = b - aMinOne;
	S = a;
	
	return P * pow( acrossAB, 3.0f ) + Q * pow( acrossAB, 2.0f ) + R * acrossAB + S;
}

#include "DarkGDK.h"
#include "Windows.h"
#include "Winuser.h"
#include <string>
#include "Math.h"

// Get desktop res
int displayX = GetSystemMetrics(SM_CXSCREEN);
int displayY = GetSystemMetrics(SM_CYSCREEN);

// defines
#define startingRes  512
#define totalOctaves 8

//#define INTERPOLATION_LINEAR
//#define INTERPOLATION_COSINE // 561
#define INTERPOLATION_CUBIC // 1459 - 1155 - 845 - 805 - 653 - 506 - 427

// Globals
float minHeight = 0.0f;
float maxHeight = 0.0f;

// Arrays
float* noiseMap[totalOctaves+1];
char   globalString[256];

// Prototypes
float MyRand( int x );
float Sample( int octaveID, int x, int y, int octaveRes );
float LinearInterpolate( float a, float b, float x );
float CubicInterpolate( float aMinOne, float a, float b, float bAddOne, float acrossAB );
float CosineInterpolate( float a, float b, float x );

void DarkGDK ( void )
{
	// Set display
	dbSetWindowPosition( 0, 0 );
	dbSetDisplayMode( displayX, displayY, 32 );
	dbSetWindowLayout( 0, 0, 0 );
	
	// Init
	dbSyncOn();
	
	// BENCHMARK ///////
	int benchmark1Start = dbTimer();
	////////////////////
	
	// Calculate weight/opacity/strength of all layers
	float strength[totalOctaves];
	float currentStrength	= 1.0;
	float maxPossibleRange	= 0.0;
	for ( int octaveID = totalOctaves; octaveID >= 1; octaveID -- )
	{
		strength[octaveID] = currentStrength;
		maxPossibleRange  += currentStrength;
		currentStrength   *= 0.75f;
	}
	
	// This method is faster than evlauating the actual min and max heights used
	minHeight = -maxPossibleRange;
	maxHeight =  maxPossibleRange;
	
	// Generate noise for all octaves
	for( int octaveID = 1; octaveID <= totalOctaves; octaveID ++ )
	{
		int octaveRes		= startingRes / pow( 2.0f, (float) (octaveID - 1) );
		noiseMap[octaveID] = new float[octaveRes * octaveRes];
		int resLoopMax = octaveRes - 1;
		//float strength = 1.0f / (float) (1 << (totalOctaves - octaveID));
		
		for( int x = 0; x <= resLoopMax; x ++ )
		{
			for( int y = 0; y <= resLoopMax; y ++ )
			{
				//noiseMap[octaveID][x+y*octaveRes] = ( ((float) dbRnd(2000)) - 1000.0f ) * 0.001;
				noiseMap[octaveID][x+y*octaveRes] = MyRand(x+y*octaveRes) * strength[octaveID];
			}
		}
	}
	
	// BENCHMARK ///////
	int benchmark1Elapsed = dbTimer() - benchmark1Start;
	int benchmark2Start = dbTimer();
	////////////////////
	
	int memblockID  = 1;
	
	dbMakeMemblock( memblockID, 12 + startingRes * startingRes * 4 );
	DWORD* memblockPtr = (DWORD*) dbGetMemblockPtr(memblockID);
	
	// Write memblock header
	memblockPtr[0] = (DWORD) startingRes;
	memblockPtr[1] = (DWORD) startingRes;
	memblockPtr[2] = (DWORD) 32;
	
	float minInv			= -minHeight;
	float heightRange		=  maxHeight - minHeight;
	float _255DivRange		=  255.0 / heightRange;
	float sResReciprocal	=    1.0 / (startingRes);
	DWORD BGRA;
	
	int resLoopMax = startingRes - 1;
	
	//sprintf( globalString, "This: %f and this %f", minHeight, maxHeight );
	//MessageBox( NULL, globalString, "Hallo", MB_OK );
	float mHeight;
	
	for( int x = 0; x <= resLoopMax; x ++ )
	{
		for( int y = 0; y <= resLoopMax; y ++ )
		{
			float fX		= ((float) x) * sResReciprocal;
			float fY		= ((float) y) * sResReciprocal;
			float height	= 0.0;
			
			for ( int octaveID = totalOctaves; octaveID >= 2; octaveID -- )
			{
				int   octaveRes	= startingRes / (1 << (octaveID - 1));
				float octavefX	= fX * (float) octaveRes;
				float octavefY	= fY * (float) octaveRes;
				int	  octaveX	= (int) octavefX;
				int	  octaveY	= (int) octavefY;
				//octavefX = fmod( octavefX, 1.0f);
				//octavefY = fmod( octavefY, 1.0f);
				octavefX = 0.0000152587890625f * (float)(((int)(octavefX * 65536.0f)) % 65536);
				octavefY = 0.0000152587890625f * (float)(((int)(octavefY * 65536.0f)) % 65536);
				
				// BICUBIC ///////////
				#ifdef INTERPOLATION_CUBIC
					float sample1  = Sample( octaveID, octaveX - 1, octaveY - 1, octaveRes );
					float sample2  = Sample( octaveID, octaveX    , octaveY - 1, octaveRes );
					float sample3  = Sample( octaveID, octaveX + 1, octaveY - 1, octaveRes );
					float sample4  = Sample( octaveID, octaveX + 2, octaveY - 1, octaveRes );
					float mHeight1 = CubicInterpolate( sample1, sample2, sample3, sample4, octavefX );

						  sample1  = Sample( octaveID, octaveX - 1, octaveY    , octaveRes );
						  sample2  = Sample( octaveID, octaveX    , octaveY    , octaveRes );
						  sample3  = Sample( octaveID, octaveX + 1, octaveY    , octaveRes );
						  sample4  = Sample( octaveID, octaveX + 2, octaveY    , octaveRes );
					float mHeight2 = CubicInterpolate( sample1, sample2, sample3, sample4, octavefX );

						  sample1  = Sample( octaveID, octaveX - 1, octaveY + 1, octaveRes );
						  sample2  = Sample( octaveID, octaveX    , octaveY + 1, octaveRes );
						  sample3  = Sample( octaveID, octaveX + 1, octaveY + 1, octaveRes );
						  sample4  = Sample( octaveID, octaveX + 2, octaveY + 1, octaveRes );
					float mHeight3 = CubicInterpolate( sample1, sample2, sample3, sample4, octavefX );

						  sample1  = Sample( octaveID, octaveX - 1, octaveY + 2, octaveRes );
						  sample2  = Sample( octaveID, octaveX    , octaveY + 2, octaveRes );
						  sample3  = Sample( octaveID, octaveX + 1, octaveY + 2, octaveRes );
						  sample4  = Sample( octaveID, octaveX + 2, octaveY + 2, octaveRes );
					float mHeight4 = CubicInterpolate( sample1, sample2, sample3, sample4, octavefX );

					mHeight  = CubicInterpolate( mHeight1, mHeight2, mHeight3, mHeight4, octavefY );
				#endif
				//////////////////////
				// COSINE ////////////
				#ifdef INTERPOLATION_COSINE
					float sample1  = Sample( octaveID, octaveX    , octaveY, octaveRes );
					float sample2  = Sample( octaveID, octaveX + 1, octaveY, octaveRes );
					float mHeightX = CosineInterpolate( sample1, sample2, octavefX );
					
						  sample1  = Sample( octaveID, octaveX    , octaveY + 1, octaveRes);
						  sample2  = Sample( octaveID, octaveX + 1, octaveY + 1, octaveRes);
					float mHeightY = CosineInterpolate( sample1, sample2, octavefX );
					
					mHeight = CosineInterpolate( mHeightX, mHeightY, octavefY );
				#endif
				/////////////////////
				// BILINEAR /////////
				#ifdef INTERPOLATION_LINEAR
					float sample1  = Sample( octaveID, octaveX    , octaveY, octaveRes );
					float sample2  = Sample( octaveID, octaveX + 1, octaveY, octaveRes );
					float mHeightX = LinearInterpolate( sample1, sample2, octavefX );
					
						  sample1  = Sample( octaveID, octaveX    , octaveY + 1, octaveRes);
						  sample2  = Sample( octaveID, octaveX + 1, octaveY + 1, octaveRes);
					float mHeightY = LinearInterpolate( sample1, sample2, octavefX );
					
					mHeight = LinearInterpolate( mHeightX, mHeightY, octavefY );
				#endif
				/////////////////////
				
				height   = height + mHeight;
			}
			// Highest frequency noise
			height   = height + Sample( 1, x, y, startingRes);
			
			height   = (height + minInv) * _255DivRange;
			int iHeight = (int) height;
			BGRA = iHeight | (iHeight<<8) | (iHeight<<16) | (255<<24);
			memblockPtr[3 + x + y * startingRes] = BGRA;
		}
	}
	
	// BENCHMARK ///////
	int benchmark2Elapsed = dbTimer() - benchmark1Start;
	////////////////////
	
	// Make image from memblock
	dbMakeImageFromMemblock( 1, memblockID );
	
	// Output image and info
	while ( LoopGDK() )
	{
		sprintf( globalString, "Time to generate noise: %ims", benchmark1Elapsed );
		dbText( 5, 5, globalString );
		sprintf( globalString, "Time to add octaves: %ims", benchmark2Elapsed );
		dbText( 5, 15, globalString );
		
		dbPasteImage( 1, 64, 64 );
		
		dbSync();
	}
	
	return;
}

float MyRand( int x )
{	 
	x = (x<<13) ^ x;
	
	return ( 1.0 - ( (x * (x * x * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0f);    
}
float Sample( int octaveID, int x, int y, int octaveRes )
{
	// Wrap
	if (x < 0)
		x += octaveRes;
	if (y < 0)
		y += octaveRes;
	
	if (x > octaveRes-1)
		x -= octaveRes;
	if (y > octaveRes-1)
		y -= octaveRes;
	
	return noiseMap[octaveID][x+y*octaveRes];
}

float LinearInterpolate( float a, float b, float x )
{
	return a * (1.0 - x) + b * x;
}
float CosineInterpolate( float a, float b, float x )
{
    x = (1.0 - cos(x * 3.141592654)) * 0.5;
    
	return a * (1.0 - x) + b * x;
}
float CubicInterpolate( float aMinOne, float a, float b, float bAddOne, float acrossAB )
{
	float T;
	
	bAddOne = (bAddOne - b) - (aMinOne - a);
	T = acrossAB * acrossAB;
	
	return bAddOne * (T * acrossAB) + ((aMinOne - a) - bAddOne) * T + (b - aMinOne) * acrossAB + a;
}

` Green Gandalf's Perlin style noise function - Version 4
` Created 25 June 2007, modified 19 February 2008.

`   Uses suggestions from following website
`      http://www.mandelbrot-dazibao.com/Perlin/Perlin1.htm
`   and especially
`      http://mrl.nyu.edu/~perlin/ (and browse to his 2002 paper)

` This version restricts the gradients to +/- 1 as suggested
` by Ken Perlin (in fact there's no need for the rnd() function at all
` - could probably just use a simple method of haphazardly swapping the signs
` of the gradients).

` All four methods work well - but simple linear filtering shows slight seams.

sync on: sync rate 0: sync
set display mode 800, 600, 32

time = timer()
randomize 140549 ` arbitrary fixed number for reproducibility

autocam off
position camera 0, 50, -300
point camera 0, 0, 0

create bitmap 1, 512, 512

global nOct = 8
global nGrid
global twoPiInv as float
twoPiInv = 0.159154943

nGrid = 2^nOct

dim weight#(nOct)

dim rawNoise(511, 511) as float

p# = 2.0
for i=0 to nOct
  weight#(i) = 1.0/p#^i
next i

dim a(nGrid, nGrid, nOct) as float
dim b(nGrid, nGrid, nOct) as float

for oct = 0 to nOct
  tGrid = 2^oct
  for i = 0 to tGrid
    for j = 0 to tGrid
      a(i, j, oct) = (rnd(1) - 0.5)* 2.0  ` random value +/- 1
      b(i, j, oct) = (rnd(1) - 0.5)* 2.0  ` random value +/- 1
                                             ` not seamless yet
    next j
  next i
next oct

` calculate raw noise and find scale factors for
` rescaling noise values to byte range
minNoise# = 10000  ` arbitrary large number out of range
maxNoise# = -10000 ` arbitrary small number out of range
for x = 0 to 511
  for y = 0 to 511
    u# = x * 0.001953125 ` value in range 0 to 1
    v# = y * 0.001953125
    rawNoise(x, y) = noise(u#, v#)
    if minNoise# > rawNoise(x, y)
      minNoise# = rawNoise(x, y)
    else
      if maxNoise# < rawNoise(x, y) then maxNoise# = rawNoise(x, y)
    endif
  next y
next x

f# = 255.0/(maxNoise# - minNoise#)

lock pixels
  for x = 0 to 511
    for y = 0 to 511
      ` convert raw noise to byte range 0 - 255
      c = (rawNoise(x, y) - minNoise#) * f#
      ` just in case
      if c<0 then c=0
      if c>255 then c=255
      dot x, y, rgb(0, c, 0) ` a nice green colour :)
    next y
  next x
unlock pixels

copy bitmap 1, 0, 0, 511, 511,  0, 144, 44, 655, 555

set current bitmap 0

time = timer() - time
repeat
  text 20, 20, "All done in "+str$(time)+" millisecs"
  sync
until spacekey()

set current bitmap 1
get image 1, 0, 0, 512, 512
save image "test v3.png", 1
end

function noiseBase(oct, x as float, y as float)
  ` x and y both in range 0 to 1
  ` find correct "tile" to use
  tGrid = 2^oct
  x# = x * tGrid
  y# = y * tGrid
  i = floor(x#)
  if i >= tGrid then i = tGrid - 1 ` just in case
  j = floor(y#)
  if j >= tGrid then j = tGrid - 1 ` just in case
  x# = x# - i  ` number in range 0 to 1
  y# = y# - j
  ` calculate interpolation weights
  ` revised Perlin weights
  remstart
  ix1# = x# * x# * x# * (6.0 * x# * x# - 15.0 * x# + 10.0)
  ix0# = 1.0 - ix1#
  iy1# = y# * y# * y# * (6.0 * y# * y# - 15.0 * y# + 10.0)
  iy0# = 1.0 - iy1#
  remend
  ` old Perlin weights - image not right at all WHY?
  remstart
  ix1# = x# * x# * (3.0 - 2.0 * x#)
  ix0# = 1.0 - ix1#
  iy1# = y# * y# * (3.0 - 2.0 * y#)
  iy0# = 1.0 - iy1#
  remend
  ` simple bilinear weights - should show slight seams
  ` because derivatives not matched correctly
  remstart
  ix0# = 1.0 - x#
  ix1# = x#
  iy0# = 1.0 - y#
  iy1# = y#
  remend
  ` simple sine weights (a variation of the Perlin idea)
`  remstart
  ix1# = x# - sin(360.0 * x#) * twoPiInv
  ix0# = 1.0 - ix1#
  iy1# = y# - sin(360.0 * y#) * twoPiInv
  iy0# = 1.0 - iy1#
`  remend
  itp0# = ix0# * iy0#
  itp1# = ix0# * iy1#
  itp2# = ix1# * iy0#
  itp3# = ix1# * iy1#
  ` calculate tangents at four corners of tile
  t0# = a(i, j, oct) * x# + b(i, j, oct) * y#
  t1# = a(i, j + 1, oct) * x# + b(i, j + 1, oct) * (y# - 1.0)
  t2# = a(i + 1 , j, oct) * (x# - 1.0) + b(i + 1, j, oct) * y#
  t3# = a(i + 1, j + 1, oct) * (x# - 1.0) + b(i + 1, j + 1, oct) * (y# - 1.0)
  ` calculate interpolated noise value
  result# = itp0#*t0# + itp1#*t1# + itp2#*t2# + itp3#*t3#
endfunction result#

function noise(x as float, y as float)
  result# = 0.0
  for oct = 0 to 8
    result# = noiseBase(oct, x, y) * weight#(oct) + result#
  next oct
endfunction result#

Load DLL "user32.dll"   , 1
displayX = Call DLL( 1 , "GetSystemMetrics" , 0 )
displayY = Call DLL( 1 , "GetSystemMetrics" , 1 )
Delete DLL 1

Set Display Mode displayX, displayY, 32, 1
Sync On

// Seed
Randomize 1338

// You can manually specify the resolution for higher detail noise
startingRes     = 512
totalOctaves    = 8

`totalOctaves    = 8
`startingRes     = 2 ^ totalOctaves

// BiCosine --------- 1946ms, 1658ms
// BiCubic  --------- 6943ms, 4967ms, 4922ms, 4903ms, 4823ms, 4672ms, 4078ms

benchmark1Start = Timer()

// Generate noise for all octaves
for octaveID = 1 to totalOctaves
    memblockID  = octaveID
    octaveRes   = startingRes / 2 ^ (octaveID - 1)
    Make Memblock memblockID, octaveRes * octaveRes * 4 // A single float
    for x = 0 to octaveRes - 1
        for y = 0 to octaveRes - 1
            Write Memblock Float memblockID, ( x + y * octaveRes ) * 4, ( Rnd(2000) - 1000 ) * 0.001
        next
    next
next

benchmark1elapsed   = Timer() - benchmark1Start

benchmark2Start = Timer()

// add all layers together
memblockID  = totalOctaves + 1

Make Memblock memblockID, 12 + startingRes * startingRes * 4

Write Memblock DWord memblockID, 0, startingRes
Write Memblock DWord memblockID, 4, startingRes
Write Memblock DWord memblockID, 8, 32

Dim strength(totalOctaves) as float
Dim res(totalOctaves)

strength# = 1.0
for octaveID = totalOctaves to 1 step -1
    strength# = strength# * 0.75
    strength(octaveID) = strength#
    res(octaveID)      = startingRes / 2 ^ (octaveID - 1)
    maxStrength#       = maxStrength# + strength#
next

minHeight#      = -maxStrength#
maxHeight#      =  maxStrength#
minInv#         = -minHeight#
heightRange#    =  maxHeight# - minHeight#
_255DivRange#   =  255.0 / heightRange#
sResReciprocal# =    1.0 / (startingRes)

startingResSubOne = startingRes - 1
for x = 0 to startingResSubOne
    for y = 0 to startingResSubOne
        x#          = x * sResReciprocal#
        y#          = y * sResReciprocal#
        height#     = 0.0
        Strength#   = 1.0

        for OctaveID = 2 to totalOctaves
            octaveRes = res(OctaveID)
            // Store the float pixel we are over on the current octave
            memblockX# = x# * octaveRes
            memblockY# = y# * octaveRes
            // Store the int version for pixel sampling
            memblockX  = memblockX#
            memblockY  = memblockY#
            // Get the local offset
            memblockX# = memblockX# mod 1.0
            memblockY# = memblockY# mod 1.0
            
            mX1 = memblockX - 1
            mX2 = memblockX
            mX3 = memblockX + 1
            mX4 = memblockX + 2
            
            mY1 = memblockY - 1
            mY2 = memblockY
            mY3 = memblockY + 1
            mY4 = memblockY + 2
            // BICUBIC ///////////
            REMSTART
            sample1#  = Sample( OctaveID, mX1, mY1, octaveRes )
            sample2#  = Sample( OctaveID, mX2, mY1, octaveRes )
            sample3#  = Sample( OctaveID, mX3, mY1, octaveRes )
            sample4#  = Sample( OctaveID, mX4, mY1, octaveRes )
            mHeight1# = CubicInterpolate( sample1#, sample2#, sample3#, sample4#, memblockX# )

            sample1#  = Sample( OctaveID, mX1, mY2, octaveRes )
            sample2#  = Sample( OctaveID, mX2, mY2, octaveRes )
            sample3#  = Sample( OctaveID, mX3, mY2, octaveRes )
            sample4#  = Sample( OctaveID, mX4, mY2, octaveRes )
            mHeight2# = CubicInterpolate( sample1#, sample2#, sample3#, sample4#, memblockX# )

            sample1#  = Sample( OctaveID, mX1, mY3, octaveRes )
            sample2#  = Sample( OctaveID, mX2, mY3, octaveRes )
            sample3#  = Sample( OctaveID, mX3, mY3, octaveRes )
            sample4#  = Sample( OctaveID, mX4, mY3, octaveRes )
            mHeight3# = CubicInterpolate( sample1#, sample2#, sample3#, sample4#, memblockX# )

            sample1#  = Sample( OctaveID, mX1, mY4, octaveRes )
            sample2#  = Sample( OctaveID, mX2, mY4, octaveRes )
            sample3#  = Sample( OctaveID, mX3, mY4, octaveRes )
            sample4#  = Sample( OctaveID, mX4, mY4, octaveRes )
            mHeight4# = CubicInterpolate( sample1#, sample2#, sample3#, sample4#, memblockX# )

            mHeight#  = CubicInterpolate( mHeight1#, mHeight2#, mHeight3#, mHeight4#, memblockY# )
            REMEND
            //////////////////////
            
            // COSINE ////////////
            REMSTART
            sample1#  = Sample( OctaveID, mX2, mY2, octaveRes )
            sample2#  = Sample( OctaveID, mX3, mY2, octaveRes )
            mHeightX# = CosineInterpolate( sample1#, sample2#, memblockX# )

            sample1#  = Sample( OctaveID, mX2, mY3, octaveRes)
            sample2#  = Sample( OctaveID, mX3, mY3, octaveRes)
            mHeightY# = CosineInterpolate( sample1#, sample2#, memblockX# )

            mHeight#  = CosineInterpolate( mHeightX#, mHeightY#, memblockY# )
            REMEND
            `REMSTART
            // Inline version
            memblockX# = (1.0 - Cos(memblockX# * 180.0)) * 0.5
            memblockY# = (1.0 - Cos(memblockY# * 180.0)) * 0.5
            
            sample1#  = Sample( OctaveID, mX2, mY2, octaveRes )
            sample2#  = Sample( OctaveID, mX3, mY2, octaveRes )
            mHeightX# = sample1# * (1.0 - memblockX#) + sample2# * memblockX#

            sample1#  = Sample( OctaveID, mX2, mY3, octaveRes)
            sample2#  = Sample( OctaveID, mX3, mY3, octaveRes)
            mHeightY# = sample1# * (1.0 - memblockX#) + sample2# * memblockX#

            mHeight#  = mHeightX# * (1.0 - memblockY#) + mHeightY# * memblockY#
            `REMEND
            //////////////////////
            
            // BILINEAR //////////
            REMSTART
            sample1#  = Sample( OctaveID, memblockX    , memblockY, octaveRes )
            sample2#  = Sample( OctaveID, memblockX + 1, memblockY, octaveRes )
            mHeightX# = LinearInterpolate( sample1#, sample2#, memblockX# )

            sample1#  = Sample( OctaveID, memblockX    , memblockY + 1, octaveRes)
            sample2#  = Sample( OctaveID, memblockX + 1, memblockY + 1, octaveRes)
            mHeightY# = LinearInterpolate( sample1#, sample2#, memblockX# )

            mHeight#  = LinearInterpolate( mHeightX#, mHeightY#, memblockY# )
            REMEND
            //////////////////////

            height#   = height# + mHeight# * strength(octaveID)
        next
        // Add octave 1
        height# = height# + Sample( 1, X, Y, startingRes ) * strength(1)
        
        height  = (height# + minInv#) * _255DivRange#
        Write Memblock DWord memblockID, 12 + ( x + y * startingRes ) * 4, height || (height<<8) || (height<<16) || (255<<24)
    next
next

benchmark2elapsed   = Timer() - benchmark2Start

// Create image
Make Image From Memblock 1, memblockID

// Loop to output result
do
    CLS

    Paste Image 1, 64, 64

    Text 5,  5, "Time to generate noise: " + Str$(benchmark1elapsed) + "ms"
    Text 5, 15, "Time to add octaves: "    + Str$(benchmark2elapsed) + "ms"

    Sync
loop

END

Function Sample( memblockID, x, y, memblockRes )

    // Wrap
    if x < 0
        x = x + memblockRes
        if y < 0             then y = y + memblockRes
        if y > memblockRes-1 then y = y - memblockRes
    else
        if x > memblockRes-1 then x = x - memblockRes
        
        if y < 0             then y = y + memblockRes
        if y > memblockRes-1 then y = y - memblockRes
    endif

    returnValue# = Memblock Float( memblockID, ( x + y * memblockRes ) * 4 )

Endfunction returnValue#

Function LinearInterpolate( a as float, b as float, x as float )

    returnValue# =  a * (1.0 - x) + b * x

Endfunction returnValue#

function CosineInterpolate( a as float, b as float, x as float )
    
    x = (1.0 - Cos(x * 180.0)) * 0.5
    
    returnValue# = a * (1.0 - x) + b * x
    
Endfunction returnValue#

Function CubicInterpolate( aMinOne as float, a as float, b as float, bAddOne as float, acrossAB as float )

    T as float
    
    bAddOne = (bAddOne - b) - (aMinOne - a)
    T = acrossAB * acrossAB
    
    returnValue# =  bAddOne * T * acrossAB + ((aMinOne - a) - bAddOne) * T + (b - aMinOne) * acrossAB + a

EndFunction returnValue#

Load DLL "user32.dll"   , 1
displayX = Call DLL( 1 , "GetSystemMetrics" , 0 )
displayY = Call DLL( 1 , "GetSystemMetrics" , 1 )
Delete DLL 1

Set Display Mode displayX, displayY, 32, 1
Sync On

// Seed
Randomize 1338

// You can manually specify the resolution for higher detail noise
startingRes     = 512
totalOctaves    = 9

`totalOctaves    = 8
`startingRes     = 2 ^ totalOctaves

// BiCosine --------- 1946ms, 1658ms, 1601ms - 9x 1783ms
// BiCubic  --------- 6943ms, 4967ms, 4922ms, 4903ms, 4823ms, 4672ms, 4078ms

// BENCHMARK //////////////
benchmark1Start = Timer()
///////////////////////////

Dim strength(totalOctaves)  as float
Dim res(totalOctaves)       as integer
Dim spacing(totalOctaves)   as integer

strength# = 1.0
for octaveID = totalOctaves to 1 step -1
    strength#          = strength# * 0.75
    strength(octaveID) = strength#
    res(octaveID)      = startingRes / 2 ^ (octaveID - 1)
    spacing(octaveID)  = startingRes / res(octaveID)
    maxStrength#       = maxStrength# + strength#
next

// Generate noise for all octaves
for octaveID = 1 to totalOctaves
    memblockID  = octaveID
    octaveRes   = startingRes / 2 ^ (octaveID - 1)
    pos         = 0
    Make Memblock memblockID, octaveRes * octaveRes * 4 // A single float
    for y = 0 to octaveRes - 1
        for x = 0 to octaveRes - 1
            Write Memblock Float memblockID, pos, ( ( Rnd(2000) - 1000 ) * 0.001 ) * strength(octaveID)
            Inc pos, 4
        next
    next
next

// BENCHMARK //////////////
benchmark1elapsed   = Timer() - benchmark1Start
benchmark2Start     = Timer()
///////////////////////////

// add all layers together
memblockID  = totalOctaves + 1

// Make memblock and write header
Make Memblock memblockID, 12 + startingRes * startingRes * 4

Write Memblock DWord memblockID, 0, startingRes
Write Memblock DWord memblockID, 4, startingRes
Write Memblock DWord memblockID, 8, 32

// Compute values
minHeight#      = -maxStrength#
maxHeight#      =  maxStrength#
minInv#         = -minHeight#
heightRange#    =  maxHeight# - minHeight#
_255DivRange#   =  255.0 / heightRange#
sResReciprocal# =    1.0 / (startingRes)

startingResSubOne = startingRes - 1
for y = 0 to startingResSubOne
    for x = 0 to startingResSubOne
        x#          = x * sResReciprocal#
        y#          = y * sResReciprocal#
        height#     = 0.0
        Strength#   = 1.0

        for OctaveID = 2 to totalOctaves
            `if (x mod spacing(octaveID))+(y mod spacing(octaveID)) = 0
            `    height# = height# + Memblock Float( OctaveID, ( x# * res(OctaveID) + y# * res(OctaveID) * res(OctaveID) ) * 4 )
            `else
            
            octaveRes = res(OctaveID)
            // Store the float pixel we are over on the current octave
            memblockX# = x# * octaveRes
            memblockY# = y# * octaveRes
            // Store the int version for pixel sampling
            memblockX  = memblockX#
            memblockY  = memblockY#
            // Get the local offset
            memblockX# = memblockX# mod 1.0
            memblockY# = memblockY# mod 1.0
            
            `if (memblockX# = 0.0 or memblockX# = 1.0) and (memblockY# = 0.0 or memblockY# = 1.0)
            `    height# = height# + Memblock Float( OctaveID, ( memblockX + memblockY * octaveRes ) * 4 )
            `else
                // BICOSINE //////////
                mX2 = memblockX
                mX3 = memblockX + 1
                mY2 = memblockY
                mY3 = memblockY + 1
                
                memblockX# = (1.0 - Cos(memblockX# * 180.0)) * 0.5
                memblockY# = (1.0 - Cos(memblockY# * 180.0)) * 0.5
                
                sample1#  = Sample( OctaveID, mX2, mY2  , octaveRes )
                sample2#  = Sample( OctaveID, mX3, mY2  , octaveRes )
                mHeightX# = sample1# * (1.0 - memblockX#) + sample2# * memblockX#
    
                sample1#  = Sample( OctaveID, mX2, mY3, octaveRes)
                sample2#  = Sample( OctaveID, mX3, mY3, octaveRes)
                mHeightY# = sample1# * (1.0 - memblockX#) + sample2# * memblockX#
    
                height# = height# + mHeightX# * (1.0 - memblockY#) + mHeightY# * memblockY#
                //////////////////////
            `endif
        next
        // Add octave 1
        height# = height# + Memblock Float( 1, ( X + Y * startingRes ) * 4 )
        
        height  = (height# + minInv#) * _255DivRange#
        Write Memblock DWord memblockID, 12 + ( x + y * startingRes ) * 4, height || (height<<8) || (height<<16) || (255<<24)
    next
next

// BENCHMARK //////////////
benchmark2elapsed   = Timer() - benchmark2Start
///////////////////////////

// Create image
Make Image From Memblock 1, memblockID

// Loop to output result
do
    CLS

    Paste Image 1, 64, 64

    Text 5,  5, "Time to generate noise: " + Str$(benchmark1elapsed) + "ms"
    Text 5, 15, "Time to add octaves: "    + Str$(benchmark2elapsed) + "ms"

    Sync
loop

END

Function Sample( memblockID, x, y, memblockRes )

    // Wrap
    if x < 0
        x = x + memblockRes
        if y < 0             then y = y + memblockRes
        if y > memblockRes-1 then y = y - memblockRes
    else
        if x > memblockRes-1 then x = x - memblockRes
        
        if y < 0             then y = y + memblockRes
        if y > memblockRes-1 then y = y - memblockRes
    endif

    returnValue# = Memblock Float( memblockID, ( x + y * memblockRes ) * 4 )

Endfunction returnValue#

Function LinearInterpolate( a as float, b as float, x as float )

    returnValue# =  a * (1.0 - x) + b * x

Endfunction returnValue#

function CosineInterpolate( a as float, b as float, x as float )
    
    x = (1.0 - Cos(x * 180.0)) * 0.5
    
    returnValue# = a * (1.0 - x) + b * x
    
Endfunction returnValue#

Function CubicInterpolate( aMinOne as float, a as float, b as float, bAddOne as float, acrossAB as float )

    T as float
    
    bAddOne = (bAddOne - b) - (aMinOne - a)
    T = acrossAB * acrossAB
    
    returnValue# =  bAddOne * T * acrossAB + ((aMinOne - a) - bAddOne) * T + (b - aMinOne) * acrossAB + a

EndFunction returnValue#

mov myVar, eax
mov eax, myVar

myVar = 1 * 2 * 3 * 4 * 5 * 6

myVar = 1 * 2 * 3 * 4 * 5 * 6

#include <math.h>
float gradientTab[TABSIZE*3];
void gradientTabInit(int seed){
 float *table = gradientTab;
 float z, r, theta;
 int i;
 srandom(seed);
 for (i=0; i < TABSIZE; i++){
  z = 1. - 2. * RANDNBR;
  /* r is radius of x,y circle */
  r = sqrtf(1 - z*z);
  /* theta is angle in (x,y) */
  theta = 2 * M_PI * RANDNBR;
  *table++ = r * cosf(theta);
  *table++ = r * sinf(theta);
  *table++ = z;
 }
}
float glattice(int ix,int iy, int iz, float fx, float fy, float fz){
 float *g = &gradientTab[INDEX(ix,iy,iz)*3];
 return g[0]*fx + g[1]*fy + g[2]*fz;
}

double Noise1D(int x){
 n = (x<<13)^x;
 return (1. -((n * (n * n * 15731 + 789221)+ 1376312589) & 7fffffff) / 1073741824.);
}
double Noise2D(int x, int y){
 n = x + y * 57;
 n = (n<<13)^n;
 return (1. -((n * (n * n * 15731 + 789221)+ 1376312589) & 7fffffff) / 1073741824.);
}
double Noise3D(int x, int y, int z){
 n = x + y * 57 + z * 131;
 n = (n<<13)^n;
 return (1. -((n * (n * n * 15731 + 789221)+ 1376312589) & 7fffffff) / 1073741824.);
}

`*** Noise functions in DBP
Function Noise1D(X as integer)
   n = (X << 13) ^X `<- I think this part isnt working correctly possibly the ^ command.
   Result# = (1.0 -((n*(n*n*15731 + 789221) + 1376312589) && 0x7fffffff) / 1073741824.0 )
EndFunction Result#
Function Noise2D(X as integer ,Y as integer)
   n = X + Y * 57
   n = (n << 13) ^n
   Result# = (1.0 -((n*(n*n*15731 + 789221) + 1376312589) && 0x7fffffff) / 1073741824.0 )
EndFunction Result#
Function Noise3D(X as integer ,Y as integer,Z as Integer)
   n = X + Y * 57 + Z * 131
   n = (n << 13) ^n
   Result# = (1.0 -((n*(n*n*15731 + 789221) + 1376312589) && 0x7fffffff) / 1073741824.0 )
EndFunction Result#

` Green Gandalf's Perlin style noise function - Version 6
` Created 25 June 2007, modified 21 February 2008.

`   Uses suggestions from following website
`      http://www.mandelbrot-dazibao.com/Perlin/Perlin1.htm
`   and especially
`      http://mrl.nyu.edu/~perlin/ (and browse to his 2002 paper)

` This version restricts the gradients to +/- 1 as suggested
` by Ken Perlin (in fact there's no need for the rnd() function at all
` - could probably just use a simple method of haphazardly swapping the signs
` of the gradients).

` All four methods work well - but simple linear filtering shows slight seams on final image.

` Attempt at an another optimized version.

sync on: sync rate 0: sync
set display mode 800, 600, 32

time0 = timer()

randomize 140549 ` arbitrary fixed number for reproducibility

autocam off
position camera 0, 50, -300
point camera 0, 0, 0

create bitmap 1, 512, 512

global nOct = 8
global twoPiInv as float
twoPiInv = 0.159154943

dim weight#(nOct)

dim rawNoise(511, 511) as float

dim a(512, 512, nOct) as float ` somewhat wasteful in space
dim b(512, 512, nOct) as float
dim interp(511) as float
dim offset(nOct, 512)
dim invSize(nOct) as float
dim tile(nOct)
dim invTile(nOct)

increment = 512
p# = 2.0
for oct = 0 to nOct
  for i = 0 to 512 step increment
    for j = 0 to 512 step increment
      a(i, j, oct) = rnd(1) * 2 - 1  ` random value +/- 1
      b(i, j, oct) = rnd(1) * 2 - 1  ` random value +/- 1
                                             ` not seamless yet
    next j
  next i
  for i = 0 to 512
    offset(oct, i) = i MOD increment
  next i
  weight#(oct) = 1.0/p#^oct
  invSize(oct) = 1.0/increment
  tile(oct) = increment
  invTile(oct) = 512/increment
  increment = increment/2
next oct

for i = 0 to 511
  x# = i/512.0
`  interp(i) = x#                                 ` simple bilinear weight
`  interp(i) = x# * x# * (3.0 - 2.0 * x#)         ` original Perlin weight
`  interp(i) = x# * x# * x# * (6.0 * x# * x# - 15.0 * x# + 10.0) ` revised Perlin weight
  interp(i) = x# - sin(360.0 * x#) * twoPiInv    ` sine weight
next i

time1 = timer()

` calculate raw noise and find scale factors for
` rescaling noise values to byte range
minNoise# = 10000  ` arbitrary large number out of range
maxNoise# = -10000 ` arbitrary small number out of range
for x = 0 to 511
  for y = 0 to 511
    rawNoise(x, y) = noise(x, y)
    if minNoise# > rawNoise(x, y)
      minNoise# = rawNoise(x, y)
    else
      if maxNoise# < rawNoise(x, y) then maxNoise# = rawNoise(x, y)
    endif
  next y
next x

time2 = timer()

f# = 255.0/(maxNoise# - minNoise#)

lock pixels
  for x = 0 to 511
    for y = 0 to 511
      ` convert raw noise to byte range 0 - 255
      c = (rawNoise(x, y) - minNoise#) * f#
      ` just in case
      if c<0 then c=0
      if c>255 then c=255
      dot x, y, rgb(0, c, 0) ` a nice green colour :)
    next y
  next x
unlock pixels

time3 = timer()

copy bitmap 1, 0, 0, 511, 511,  0, 144, 44, 655, 555

set current bitmap 0

time = timer()

repeat
  text 20, 20, "All done in "+str$(time-time0)+" ms"
  text 20, 40, "Table set-up = "+str$(time1-time0)+" ms"
  text 20, 60, "Raw noise set-up = "+str$(time2-time1)+" ms"
  text 20, 80, "Noise rescale and image creation = "+str$(time3-time2)+" ms"
  sync
until spacekey()

set current bitmap 1
get image 1, 0, 0, 512, 512
save image "test v3.png", 1
end

function noiseBase(oct, x, y)
  ` x and y both in range 0 to 511
  ` find x and y offsets within tile
  xd = offset(oct, x)
  yd = offset(oct, y)
  w = tile(oct)
  x# = xd * invSize(oct)  ` number in range 0 to 1
  y# = yd * invSize(oct)
  ` look-up interpolation weights
  ix# = interp(xd*invTile(oct))
  iy# = interp(yd*invTile(oct))
  ` calculate tangents at four corners of tile
  i = x - xd
  j = y - yd
  t0# = a(i, j, oct) * x# + b(i, j, oct) * y#
  t1# = a(i, j + w, oct) * x# + b(i, j + w, oct) * (y# - 1.0)
  t2# = a(i + w , j, oct) * (x# - 1.0) + b(i + w, j, oct) * y#
  t3# = a(i + w, j + w, oct) * (x# - 1.0) + b(i + w, j + w, oct) * (y# - 1.0)
  ` calculate interpolated noise value in two stages
  tA# = t0# + (t2# - t0#) * ix#
  tB# = t1# + (t3# - t1#) * ix#
  result# = tA# + (tb# - tA#) * iy#
endfunction result#

function noise(x, y)
  result# = 0.0
`  result# = noiseBase(2, x, y) ` used for testing single octave instead of next 3 lines
  for oct = 0 to nOct
    result# = noiseBase(oct, x, y) * weight#(oct) + result#
  next oct
endfunction result#

      a(i, j, oct) = rnd(1) * 2 - 1  ` random value +/- 1
      b(i, j, oct) = rnd(1) * 2 - 1  ` random value +/- 1

      a(i, j, oct) = rnd(1000) - 500 ` random value in range -500 to 500
      b(i, j, oct) = rnd(1000) - 500 ` random value in range -500 to 500

` Green Gandalf's fractal noise texture demo v2
`
` Created 13 Sep 2006, modified 22 February 2008.

` Creates a noise image using the diamond-square algorithm.
` Each colour channel contains its own noise image.
` Produces THREE noise images far faster than Perlin noise
` produces ONE.

set display mode 800,600,32
set text opaque
sync on: sync rate 60: sync
autocam off

time = timer()

set cursor 20,20

n=9: dispersion#=512: k=0 ` n=8 gives a 256x256 bitmap

maxgrid=2^n: g=maxgrid
dim col(3,maxgrid,maxgrid) as float
base as dword
global cx#: global cy#: global cz#

im=1: gosub makeImage ` prepare initial image

end

make_bitmap:
  create bitmap im,maxgrid+1,maxgrid+1
  lock pixels
    for i=0 to maxgrid
      for j=0 to maxgrid
        gosub dotij
      next j
    next i
  unlock pixels
  ` save base image
  get image 1, 0, 0, maxgrid, maxgrid
  if file exist("baseImage.png") then delete file "baseImage.png"
  save image "baseImage.png",1
  set current bitmap 0
  cls
  copy bitmap im,0
  time = timer() - time
  center text screen width()/2,screen height()/2+20,"Fractal image created in "+str$(time)+" ms"
  center text screen width()/2,screen height()/2+40,"Press any key to exit"
  sync
  wait key
return

dotij:
  dot i,j,rgb(col(1,i,j),col(2,i,j),col(3,i,j))
return

diamond_step:
  i=mid: i1=i-mid: i2=i+mid
  while i<maxgrid
    j=mid: j1=j-mid: j2=j+mid
    while j<maxgrid
      av#=(col(r,i1,j1)+col(r,i1,j2)+col(r,i2,j1)+col(r,i2,j2))/4.0
      ` calculate random values between -1 and 1
      u#=rnd(16384)/8192.0-1.0
      col(r,i,j)=av#+u#*d#`+m#
      inc j,g: inc j1,g: inc j2,g
    endwhile
    inc i,g: inc i1,g: inc i2,g
  endwhile
return

square_step:
  i=0: i1=i-mid: i2=i+mid: js=0
  while i<maxgrid
    js=mid-js ` toggle start values of j loop
    j=js: j1=j-mid: j2=j+mid
    while j<maxgrid`+1
      av#=0
      if i1<0  ` check for need to wrap around i value
        inc av#,col(r,i2,j)+col(r,i2,j)
      else
        if i2>maxgrid
          inc av#,col(r,i1,j)+col(r,i1,j)
        else
          inc av#,col(r,i1,j)+col(r,i2,j)
        endif
      endif
      if j1<0  ` check for need to wrap around j value
        inc av#,col(r,i,j2)+col(r,i,j2)
      else
        if j2>maxgrid
          inc av#,col(r,i,j1)+col(r,i,j1)
        else
          inc av#,col(r,i,j1)+col(r,i,j2)
        endif
      endif
      av#=av#/4
      ` calculate random value between -1 and 1
      u#=rnd(16384)/8192.0-1.0
      col(r,i,j)=av#+u#*d#`+m#
      col(r,maxgrid,j)=col(r,0,j) ` copy opposite edge
      inc j,g: inc j1,g: inc j2,g
    endwhile
    if j=maxgrid then col(r,i,j)=col(r,i,0) ` copy opposite edge
    inc i,mid: inc i1,mid: inc i2,mid
  endwhile
return

makeImage:
  ` Produces a fractal image based on the "diamond-square algorithm"
  ` described in following website
  `          http://www.gameprogrammer.com/fractal.html
  cls
  randomize timer()
  for r=1 to 3  ` process each colour component separately
    g=maxgrid
    d#=dispersion#
    for i=0 to maxgrid
      for j=0 to maxgrid
        col(r,i,j)=k
      next j
    next i

    ` main loop
    while g>1
      mid=g/2

      ` diamond step - calculates new diamond corners from squares
      gosub diamond_step

      ` square step - calculates new square corners from diamonds
      gosub square_step

      d#=d#/2.0: m#=m#/2.0: g=g/2
    endwhile

    ` now scale heightmap values to byte range
    min#=col(r,0,0): max#=min#
    for i=0 to maxgrid
      for j=0 to maxgrid
        u#=col(r,i,j)
        if u#<min#
          min#=u#
        else
          if u#>max# then max#=u#
        endif
      next j
    next i
    max#=256.0/(max#-min#)
    for i=0 to maxgrid
      for j=0 to maxgrid
        temp=int((col(r,i,j)-min#)*max#)
        ` range check for byte just in case
        if temp<0
          col(r,i,j)=0
        else
          if temp>255
            col(r,i,j)=255
          else
            col(r,i,j)=temp
          endif
        endif
      next j
    next i
  next r

  gosub make_bitmap

return