Improved Perlin noise

Here is an implementation of the improved Perlin noise, for four dimensions. To know more about Perlin noise, see the talk Making noise.

First, a set of gradients, aiming at the center of each edge of a hypercube, is created.

static const char gradient[32][4] =
 {
  { 1, 1, 1, 0}, { 1, 1, 0, 1}, { 1, 0, 1, 1}, { 0, 1, 1, 1},
  { 1, 1, -1, 0}, { 1, 1, 0, -1}, { 1, 0, 1, -1}, { 0, 1, 1, -1},
  { 1, -1, 1, 0}, { 1, -1, 0, 1}, { 1, 0, -1, 1}, { 0, 1, -1, 1},
  { 1, -1, -1, 0}, { 1, -1, 0, -1}, { 1, 0, -1, -1}, { 0, 1, -1, -1},
  {-1, 1, 1, 0}, {-1, 1, 0, 1}, {-1, 0, 1, 1}, { 0, -1, 1, 1},
  {-1, 1, -1, 0}, {-1, 1, 0, -1}, {-1, 0, 1, -1}, { 0, -1, 1, -1},
  {-1, -1, 1, 0}, {-1, -1, 0, 1}, {-1, 0, -1, 1}, { 0, -1, -1, 1},
  {-1, -1, -1, 0}, {-1, -1, 0, -1}, {-1, 0, -1, -1}, { 0, -1, -1, -1},
};

An array dedicated to store random numbers, and a function to fill it are also written. This function has to be called once before using Perlin noise.

static int permut[256];

//
// Function initializing the Perlin Noise
//
void InitNoise ()
{
 // Permutations
 for (unsigned int i = 0; i < 256; i++)
  permut[i] = rand () & 0xff;
}

We also need a function to return a random number for any interger coordinates, some functions to make interpolations, and a function to return the dot product of two vectors. Note that interpolation is a 5th order polynom, and the vector elements are -1, 0, or 1 (these are part of the improvement over original Perlin noise).

//
// Function finding out the gradient corresponding to the coordinates
//
static int Indice (const int i,
          const int j,
          const int k,
          const int l)
{
 return (permut[(l + permut[(k + permut[(j + permut[i & 0xff])
                     & 0xff])
               & 0xff])
         & 0xff]
     & 0x1f);
}

//
// Functions computing the dot product of the vector and the gradient
//
static float Prod (const float a, const char b)
{
 if (b > 0)
  return a;
 if (b < 0)
  return -a;
 return 0;
}

static float Dot_prod (const float x1, const char x2,
            const float y1, const char y2,
            const float z1, const char z2,
            const float t1, const char t2)
{
 return (Prod (x1, x2) + Prod (y1, y2) + Prod (z1, z2) + Prod (t1, t2));
}

static float LinearInterpolation (const float start,
                 const float end,
                 const float state)
{
 return start + state * (end - start);
}

static float Spline5 (const float state)
{
 // 3x^2 + 2x^3 is not as good as 6x^5 - 15x^4 + 10x^3
 const float square = state * state;
 const float cubic = square * state;

 return cubic * (6 * square - 15 * state + 10);
}

At last, the Perlin noise function gets the eight gradients corresponding to the given coordinates, and interpolate them to find the weight at those coordinates.

//
// Noise function, returning the Perlin Noise at a given point
//
float Noise (const float x,
       const float y,
       const float z,
       const float t)
{
 // The unit hypercube containing the point
 const int x1 = (int) (x > 0 ? x : x - 1);
 const int y1 = (int) (y > 0 ? y : y - 1);
 const int z1 = (int) (z > 0 ? z : z - 1);
 const int t1 = (int) (t > 0 ? t : t - 1);
 const int x2 = x1 + 1;
 const int y2 = y1 + 1;
 const int z2 = z1 + 1;
 const int t2 = t1 + 1;

 // The 16 corresponding gradients
 const char * g0000 = gradient[Indice (x1, y1, z1, t1)];
 const char * g0001 = gradient[Indice (x1, y1, z1, t2)];
 const char * g0010 = gradient[Indice (x1, y1, z2, t1)];
 const char * g0011 = gradient[Indice (x1, y1, z2, t2)];
 const char * g0100 = gradient[Indice (x1, y2, z1, t1)];
 const char * g0101 = gradient[Indice (x1, y2, z1, t2)];
 const char * g0110 = gradient[Indice (x1, y2, z2, t1)];
 const char * g0111 = gradient[Indice (x1, y2, z2, t2)];
 const char * g1000 = gradient[Indice (x2, y1, z1, t1)];
 const char * g1001 = gradient[Indice (x2, y1, z1, t2)];
 const char * g1010 = gradient[Indice (x2, y1, z2, t1)];
 const char * g1011 = gradient[Indice (x2, y1, z2, t2)];
 const char * g1100 = gradient[Indice (x2, y2, z1, t1)];
 const char * g1101 = gradient[Indice (x2, y2, z1, t2)];
 const char * g1110 = gradient[Indice (x2, y2, z2, t1)];
 const char * g1111 = gradient[Indice (x2, y2, z2, t2)];

 // The 16 vectors
 float dx1 = x - x1;
 float dx2 = x - x2;
 float dy1 = y - y1;
 float dy2 = y - y2;
 float dz1 = z - z1;
 float dz2 = z - z2;
 float dt1 = t - t1;
 float dt2 = t - t2;

 // The 16 dot products
 const float b0000 = Dot_prod(dx1, g0000[0], dy1, g0000[1],
                dz1, g0000[2], dt1, g0000[3]);
 const float b0001 = Dot_prod(dx1, g0001[0], dy1, g0001[1],
                dz1, g0001[2], dt2, g0001[3]);
 const float b0010 = Dot_prod(dx1, g0010[0], dy1, g0010[1],
                dz2, g0010[2], dt1, g0010[3]);
 const float b0011 = Dot_prod(dx1, g0011[0], dy1, g0011[1],
                dz2, g0011[2], dt2, g0011[3]);
 const float b0100 = Dot_prod(dx1, g0100[0], dy2, g0100[1],
                dz1, g0100[2], dt1, g0100[3]);
 const float b0101 = Dot_prod(dx1, g0101[0], dy2, g0101[1],
                dz1, g0101[2], dt2, g0101[3]);
 const float b0110 = Dot_prod(dx1, g0110[0], dy2, g0110[1],
                dz2, g0110[2], dt1, g0110[3]);
 const float b0111 = Dot_prod(dx1, g0111[0], dy2, g0111[1],
                dz2, g0111[2], dt2, g0111[3]);
 const float b1000 = Dot_prod(dx2, g1000[0], dy1, g1000[1],
                dz1, g1000[2], dt1, g1000[3]);
 const float b1001 = Dot_prod(dx2, g1001[0], dy1, g1001[1],
                dz1, g1001[2], dt2, g1001[3]);
 const float b1010 = Dot_prod(dx2, g1010[0], dy1, g1010[1],
                dz2, g1010[2], dt1, g1010[3]);
 const float b1011 = Dot_prod(dx2, g1011[0], dy1, g1011[1],
                dz2, g1011[2], dt2, g1011[3]);
 const float b1100 = Dot_prod(dx2, g1100[0], dy2, g1100[1],
                dz1, g1100[2], dt1, g1100[3]);
 const float b1101 = Dot_prod(dx2, g1101[0], dy2, g1101[1],
                dz1, g1101[2], dt2, g1101[3]);
 const float b1110 = Dot_prod(dx2, g1110[0], dy2, g1110[1],
                dz2, g1110[2], dt1, g1110[3]);
 const float b1111 = Dot_prod(dx2, g1111[0], dy2, g1111[1],
                dz2, g1111[2], dt2, g1111[3]);

 // Then the interpolations, down to the result
 dx1 = Spline5 (dx1);
 dy1 = Spline5 (dy1);
 dz1 = Spline5 (dz1);
 dt1 = Spline5 (dt1);

 const float b111 = LinearInterpolation (b1110, b1111, dt1);
 const float b110 = LinearInterpolation (b1100, b1101, dt1);
 const float b101 = LinearInterpolation (b1010, b1011, dt1);
 const float b100 = LinearInterpolation (b1000, b1001, dt1);
 const float b011 = LinearInterpolation (b0110, b0111, dt1);
 const float b010 = LinearInterpolation (b0100, b0101, dt1);
 const float b001 = LinearInterpolation (b0010, b0011, dt1);
 const float b000 = LinearInterpolation (b0000, b0001, dt1);

 const float b11 = LinearInterpolation (b110, b111, dz1);
 const float b10 = LinearInterpolation (b100, b101, dz1);
 const float b01 = LinearInterpolation (b010, b011, dz1);
 const float b00 = LinearInterpolation (b000, b001, dz1);

 const float b1 = LinearInterpolation (b10, b11, dy1);
 const float b0 = LinearInterpolation (b00, b01, dy1);

 return LinearInterpolation (b0, b1, dx1);
}

The code could probably be shortened and optimised, but at least it is clear. By the way, there is no warranty this code is bug free. :-)

Comments are welcomed of course.

This is a very good contribution. Thanks for sharing it!

You are welcome indeed. :-)

Isn't the GPL a bit radical for such a small code snippet? Anyone using it would have their entire project "infected" by the GPL, which may deter quite a lot people from using it (well, perhaps they would use it and not comply with the license). Perhaps a BSD-style license is more applicable.

Originally Posted by Kippesoep

Isn't the GPL a bit radical for such a small code snippet?

A bit, I must admit. :-)

But I think this is fair to give some free code and expect users to do so in return. By the way, it should result very easy to re-write a different (and better ^_^) implementation of this.

The problem with BSD like licences is that would allow programmers to simply copy-paste it in some commercial proprietary software, and this is just not the purpose of sharing it.

"The problem with BSD like licences is that would allow programmers to simply copy-paste it in some commercial proprietary software, and this is just not the purpose of sharing it."

The GPL license is suited for applications and libraries and not for individual code files. This code is probably around 500 lines long so I don't think it's reasonable at all to ask rest of the program in return. Also the success of the commercial product is hardly dependend on a one small code snippet.

People who have made decision to release their program under GPL are maybe ok with this, but for anyone else this code is totally useless. GPL enforces that the same license is being applied for the whole program, which is definitely an overkill for this. It would be certainly enough to have a license that requires people publish their modifications only to the code above.

I don't think the code snapshot even should allow GPL'ed entries.



With all respect,


Juhani

Originally Posted by juhnu

The GPL license is suited for applications and libraries and not for individual code files. This code is probably around 500 lines long so I don't think it's reasonable at all to ask rest of the program in return.

I agree with this argument, so I remove the statement. But I still think the purpose here is not to provide ready to paste code for employed programmers (a bit exaggerated, but you get the idea ^_^). Thank you for your points.

very nice code, thanks!

i did a quick n dirthy optimization if you want to use the code with a constant t and z during frame (you only need to change t and z between two frames to get an animated "plasma")

----

// some global vars to speed up the calculation if t and z are constant during one frame
int g_z1;
int g_t1;
int g_z2;
int g_t2;
float g_dz1;
float g_dt1;
float g_dz2;
float g_dt2;
float g_Sdz1;
float g_Sdt1;

// must be called ONCE before NemaSpeed_Noise is called in a loop
void NemaSpeed_Precalc(const float z, const float t)
{
	g_z1 = (int) (z > 0 ? z : z - 1);
	g_t1 = (int) (t > 0 ? t : t - 1);
	g_z2 = g_z1 + 1;
	g_t2 = g_t1 + 1;

	g_dz1 = z - g_z1;
	g_dz2 = z - g_z2;
	g_dt1 = t - g_t1;
	g_dt2 = t - g_t2;

	g_Sdz1 = Spline5 (g_dz1);
	g_Sdt1 = Spline5 (g_dt1);
}

// NEMA-accelerated Noise function for constant z and t during one frame, returning the Perlin Noise at a given point
float NemaSpeed_Noise (const float x, const float y)
{
	// The unit hypercube containing the point
	const int x1 = (int) (x > 0 ? x : x - 1);
	const int y1 = (int) (y > 0 ? y : y - 1);
	const int x2 = x1 + 1;
	const int y2 = y1 + 1;

	// The 16 corresponding gradients
	const char * g0000 = g_gradient[Indice (x1, y1, g_z1, g_t1)];
	const char * g0001 = g_gradient[Indice (x1, y1, g_z1, g_t2)];
	const char * g0010 = g_gradient[Indice (x1, y1, g_z2, g_t1)];
	const char * g0011 = g_gradient[Indice (x1, y1, g_z2, g_t2)];
	const char * g0100 = g_gradient[Indice (x1, y2, g_z1, g_t1)];
	const char * g0101 = g_gradient[Indice (x1, y2, g_z1, g_t2)];
	const char * g0110 = g_gradient[Indice (x1, y2, g_z2, g_t1)];
	const char * g0111 = g_gradient[Indice (x1, y2, g_z2, g_t2)];
	const char * g1000 = g_gradient[Indice (x2, y1, g_z1, g_t1)];
	const char * g1001 = g_gradient[Indice (x2, y1, g_z1, g_t2)];
	const char * g1010 = g_gradient[Indice (x2, y1, g_z2, g_t1)];
	const char * g1011 = g_gradient[Indice (x2, y1, g_z2, g_t2)];
	const char * g1100 = g_gradient[Indice (x2, y2, g_z1, g_t1)];
	const char * g1101 = g_gradient[Indice (x2, y2, g_z1, g_t2)];
	const char * g1110 = g_gradient[Indice (x2, y2, g_z2, g_t1)];
	const char * g1111 = g_gradient[Indice (x2, y2, g_z2, g_t2)];

	// The 16 vectors
	float dx1 = x - x1;
	float dx2 = x - x2;
	float dy1 = y - y1;
	float dy2 = y - y2;

	// The 16 dot products
	const float b0000 = Dot_prod(dx1, g0000[0], dy1, g0000[1], g_dz1, g0000[2], g_dt1, g0000[3]);
	const float b0001 = Dot_prod(dx1, g0001[0], dy1, g0001[1], g_dz1, g0001[2], g_dt2, g0001[3]);
	const float b0010 = Dot_prod(dx1, g0010[0], dy1, g0010[1], g_dz2, g0010[2], g_dt1, g0010[3]);
	const float b0011 = Dot_prod(dx1, g0011[0], dy1, g0011[1], g_dz2, g0011[2], g_dt2, g0011[3]);
	const float b0100 = Dot_prod(dx1, g0100[0], dy2, g0100[1], g_dz1, g0100[2], g_dt1, g0100[3]);
	const float b0101 = Dot_prod(dx1, g0101[0], dy2, g0101[1], g_dz1, g0101[2], g_dt2, g0101[3]);
	const float b0110 = Dot_prod(dx1, g0110[0], dy2, g0110[1], g_dz2, g0110[2], g_dt1, g0110[3]);
	const float b0111 = Dot_prod(dx1, g0111[0], dy2, g0111[1], g_dz2, g0111[2], g_dt2, g0111[3]);
	const float b1000 = Dot_prod(dx2, g1000[0], dy1, g1000[1], g_dz1, g1000[2], g_dt1, g1000[3]);
	const float b1001 = Dot_prod(dx2, g1001[0], dy1, g1001[1], g_dz1, g1001[2], g_dt2, g1001[3]);
	const float b1010 = Dot_prod(dx2, g1010[0], dy1, g1010[1], g_dz2, g1010[2], g_dt1, g1010[3]);
	const float b1011 = Dot_prod(dx2, g1011[0], dy1, g1011[1], g_dz2, g1011[2], g_dt2, g1011[3]);
	const float b1100 = Dot_prod(dx2, g1100[0], dy2, g1100[1], g_dz1, g1100[2], g_dt1, g1100[3]);
	const float b1101 = Dot_prod(dx2, g1101[0], dy2, g1101[1], g_dz1, g1101[2], g_dt2, g1101[3]);
	const float b1110 = Dot_prod(dx2, g1110[0], dy2, g1110[1], g_dz2, g1110[2], g_dt1, g1110[3]);
	const float b1111 = Dot_prod(dx2, g1111[0], dy2, g1111[1], g_dz2, g1111[2], g_dt2, g1111[3]);

	// Then the interpolations, down to the result
	dx1 = Spline5 (dx1);
	dy1 = Spline5 (dy1);

	const float b111 = LinearInterpolation (b1110, b1111, g_Sdt1);
	const float b110 = LinearInterpolation (b1100, b1101, g_Sdt1);
	const float b101 = LinearInterpolation (b1010, b1011, g_Sdt1);
	const float b100 = LinearInterpolation (b1000, b1001, g_Sdt1);
	const float b011 = LinearInterpolation (b0110, b0111, g_Sdt1);
	const float b010 = LinearInterpolation (b0100, b0101, g_Sdt1);
	const float b001 = LinearInterpolation (b0010, b0011, g_Sdt1);
	const float b000 = LinearInterpolation (b0000, b0001, g_Sdt1);

	const float b11 = LinearInterpolation (b110, b111, g_Sdz1);
	const float b10 = LinearInterpolation (b100, b101, g_Sdz1);
	const float b01 = LinearInterpolation (b010, b011, g_Sdz1);
	const float b00 = LinearInterpolation (b000, b001, g_Sdz1);

	const float b1 = LinearInterpolation (b10, b11, dy1);
	const float b0 = LinearInterpolation (b00, b01, dy1);

	return LinearInterpolation (b0, b1, dx1);
}

and your loop could look somehow like this:

int* pDest = NULL; // Replace NULL with your 32 bit buffer where you want to render the PerlinNoise into. Buffer must be of size m_nXRes*m_nYRes*sizeof(int);

int n = 0;
float z = 0.0f; // replace 0.0f with your animated z function
float t = 0.0f; // replace 0.0f with your animated t function

float zoom = 1.0f; // replace 1.0f with your animated zoom-function
float yStep = zoom/(float)(m_nYRes);
float xStep = zoom/(float)(m_nXRes);

int col, nX, nY;

NemaSpeed_Precalc(z,t);

float x, y = 0.0f;
for (nY = 0; nY < m_nYRes; nY++)
{
	x = 0.0f;
	for (nX = 0; nX < m_nXRes; nX++)
	{
//		col = 128+(int)(Noise(x,y,z,t)*116.0f);
		col = 128+(int)(NemaSpeed_Noise(x,y)*116.0f);
		pDest[n++] = (col << 16) | (col <<8) | col;
		x += xStep;
	}
	y += yStep;
}

have fun!

While googling for new appearances of Perlin Noise out of lazyness, I discovered this thread.

Thanks for that implementation!

And everybody who has a problem with the GPL: Just don't use this piece of code, if he wouldn't have released it at all, there wouldn't be a difference (except for that you can still study the code, wink wink). It seems that Non-GPL-People often beef the GPLers far more often than the GPL-People beef the Non-GPLers. Again: If you have a problem with the GPL, go away and don't suck on this thread. The GPL people won't annoy you so much about making your proprietary software free as you bother them.

I'm sorry, but what a load of crap. It isn't about who is bullshiting whom. It isn't even about GPL. It's about how much the code contributes to devmaster. And with a restrictive license (whether that be GPL or the PhresnelLikesToKickOldThreads license) it doesn't do that much good. If you don't want to share your code on non-intrusive terms, don't share it via devmaster.

"... or the PhresnelLikesToKickOldThreads license" - .oisyn

I hate that license!

DotOisynLovesCamelCaseAndGlue:

The thread wasn't locked. Also, I bet there are actually some people around on devmaster who like copyleft licensery, and who release stuff under such terms, so why not? I must have missed the clause about the "anything that isn't useful to at least 90% of people" embargo when I registered up to devmaster.

And as long as there isn't that embargo, wouldn't bashing on copyleft licensed releases be considered vigilante justice by a lawyer?

And even if the code is pretty short, if I read it right, isn't the category saying "Code [...] Discussion"? So let's discuss about the code, not the license.

edit: Oh, and please feel free to generalise my statements to "If you don't like ..., go away and don't bother bashing.". I wasn't actually talking specifically about GPL in that phrase, but on every license, may it be proprietary or permissive/restrictive open source. (believe it or not, I go away instead of pressuring people to release under terms that I would like.)


monjardin:
Aha.

Hmm. I am loathe to jump on a zombie thread.

The GPL is "restrictive", but I balk at calling it intrusive. It is, obviously, more intrusive that dropping snippets into the public domain. OTOH, it is less restrictive than most proprietary licenses.

However, the GPL permits you to create and distribute an "aggregate", even when the licenses of the other software are non-free or GPL-incompatible. The problem is that an "aggregate" is not defined anywhere in any technical terms and trying to ensure you don't run afowl of it is hard. See this link for a FAQ answer that addresses this.

Originally Posted by phresnel

(believe it or not, I go away instead of pressuring people to release under terms that I would like.)

Yet you are pressuring people to stop advising (I'm sorry, but I don't see any bashful post regarding the GPL in this topic) not to use the GPL for small code fragments. A sound advise, which has been taken for granted by the poster of the snippet. So, what exactly is your problem?

And the fact that a thread isn't locked doesn't mean it's ok to kick it if the post has nothing to do with the discussion at hand.