Unexpected values returned by C++ noise generation functions

algorithmic-randomnessalgorithms

I'm porting the C++ libnoise library into another language. While doing so, I came across a function which appears to return inappropriate values, at least according to the original documentation from the library. I originally thought this was a mistake on my end while porting, but I ran the original code through a C++ interpreter and get the same results as my ported code.

For reference, I have been using, and will be linking to a modern fork on GitHub. However, I've checked and made sure the functions are the same in either case, so I know its not a matter of the fork messing something up. The fork basically just modernizes the C++ code, but leaves implementation the same.

The function in question is double GradientNoise3D (double fx, double fy, double fz, int ix, int iy, int iz, int seed = 0). Here is the relevant header file and here is the implementation. The documentation claims that GradientNoise3D should be returning values between -1.0 and +1.0. However, in my testing, values range roughly between -3.3 to +3.3. This seems to trickle down to GradientCoherentNoise3D, which uses GradientNoise3D as a base function to produce noise, and similarly, claims to produce values between -1.0 and +1.0, but ends up producing values in about the same range as the base GradientNoise3D function.

It should be noted that the other functions in the noisegen namespace function as advertised. It's only the gradient noise functions which have this issue.

So now I'm trying to figure out:

  1. Is this an error in the documentation, or in the code? On one hand, a documentational error would seem more likely, in that the libnoise library is very old and seems to be trusted quite a bit. You'd think an actual code error would have been sniffed out and fixed. However, as far as my understanding goes, -3.3 to +3.3 is a very strange range for a coherent noise function, with most coherent noise functions operating within the -1.0 to +1.0 range advertised. So that would indicate a code error.
  2. Assuming this is actually not intentional behavior, what needs to be changed to provide answers within the typical -1.0 to +1.0 range? I tried changing the * 2.12 at the end of the return line in GradientNoise3D to a / 1.55... and while that more or less seems to work, feels a bit like I'm just arbitrarily brute-forcing a "solution", and was hoping to get some help from people who have more understanding of this stuff. Especially because I'm not even sure the actual theoretical limit is +/-3.3, but rather maybe some number just barely short.

Note that the g_randomVectors object is a 1024 length list of random values between -1.0 and +1.0, described in vectorable.h

Thanks to anyone who can provide some insight on this!

Best Answer

Had a look at the header and sources files you linked, those two aren't consistent with one another... Also as a disclaimer, I have no clue what a gradient-noise is supposed to be.

The supposedly computed value versus the computed value

First some notations. Let $F$ denote the point with floating point coordinates (fx, fy, fz), while $I$ is the point with integer coordinates (ix, iy, iz). According to the header file, GradientNoise3D generates a random unit vector $\vec u$, then returns the dot product: $\langle F,\,I+\vec u\rangle$

However if you look at the source, what is actually computed is: $\langle F-I,\,\vec u\rangle \times 2.12$

Range of each value

Value based on the source: Interpreting the pre-requisites as $\lvert fx-ix\rvert \le 1$, $\lvert fy-iy\rvert \le 1$, and $\lvert fz-iz\rvert \le 1$, and after checking that the $256$ possible values of $\vec u$ indeed satisfy $\left\|\vec u\right\| = 1$ (more or less), then \begin{align*} \left\lvert \langle F-I,\,\vec u\rangle \times 2.12 \right\rvert &\le 2.12\times\left\| F-I\right\|\times\left\|\vec u\right\| \\ &= 2.12\times\left\|F-I\right\|\\ &\le 2.12\sqrt 3 \\ &\approx 3.6719\ldots \end{align*} To achieve that worst case, you need to have equality in all three pre-requisites, and you'd also have to be unlucky enough to end up with $\vec u$ collinear with $F-I$, so I think it's unlikely to happen.

Value based on the header: Completely unbounded. If you pick a point $F$ that is very far from the origin, then $\langle F,\, I+\vec u\rangle$ is basically an approximation of $\left\|F\right\|^2$.

Possible correction

Again, I have no clue what a gradient noise is, but if your purpose is to get a randomized value between $-1$ and $1$, and the way $\vec u$ is computed is enough to guarantee the "random" part, you could choose to return $$ \cos(F-I,\vec u) = \frac{\langle F-I,\,\vec u\rangle}{\left\|F-I\right\|\times\left\|\vec u\right\|} $$ Since $\vec u$ already has unit norm, you could just norm $F-I$ before computing the dot product. No clue about how nice (or bad) that noise value is.

Edit: Clarification

The dot product has a nice geometric interpretation (wikipedia, mathworld). My last formula is based on that, though I guess my favourite dot product notation can throw some people off:

  • $\langle F-I,\,\vec u\rangle$ denotes the dot product between vector $F-I$ and vector $\vec u$,
  • $\left\|F-I\right\|$ denotes the norm of vector $F-I$,
  • $\left\|\vec u\right\|$ likewise denotes the norm of vector $\vec u$.

By assumption/construction, we already know that $\left\|\vec u\right\|=1$ so we can forget about that. Due to the dot product properties, and because the original c++ code already computes the three coordinates of vector $F-I$, I suggested to return $$\left\langle \frac{F-I}{\left\|F-I\right\|},\,\vec u\right\rangle$$ So all you'd have to do is norm $F-I$ before computing the dot product:

double xvPoint = (fx - (double)ix);
double yvPoint = (fy - (double)iy);
double zvPoint = (fz - (double)iz);

// modification starts below
double invNorm = 1.0 / std::sqrt( xvPoint*xvPoint + yvPoint*yxPoint + zvPoint*zvPoint );
xvPoint *= invNorm;
yvPoint *= invNorm;
zvPoint *= invNorm;

return ((xvGradient * xvPoint) + (yvGradient * yvPoint) + (zvGradient * zvPoint));
//return ((xvGradient * xvPoint) + (yvGradient * yvPoint) + (zvGradient * zvPoint)) * 2.12;
Related Question