Edit on 2023/04/13: The final version, extended to any dimension, is available in this article.

In a previous article I have introduced a noise generation algorithm based on value noise and Bezier curves. That noise suffered from directional bias, then I've modified it to solve this problem. In the present article I'm introducing a new, directionally unbiased, version of the B-Noise algorithm.

The bias in the previous version came from the choice of using two Bezier curves per axis, each using one axis as input. For a given value in one axis, one of the curve for each axis will have a fixed value whatever the value of the other axis. It creates a fixed point aligned with the axes, hence with the underlying grid, which the funtion \(F()\) fails to hide. Choosing two Bezier surfaces instead, one per mapped coordinate, these coordinates vary with both input coordinates. This decorrelates completely the inputs and the underlying grid, avoiding directional bias. Of course, the new version is still continuous, differentiable, tilable. It also stays simple to implement, and efficient to compute.

The pseudo-code becomes:

STRUCT BNoise: dim: int, number of vertices per axis (square grid) map: array of dim x dim floats in [0,1], value of each vertex in the noise grid surf: array of 2 x ((dim+1)*3+4) x ((dim+1)*3+4) floats in [0,dim-1], controls of the 2 Bezier surfaces FUNCTION PRN(): return a pseudo-random number in [0,1] FUNCTION BNoiseInit( bnoise: BNoise, the noise to initialise, order: int, the order of the noise): bnoise.dim = order + 1 nbCtrl = 3 * bnoise.dim + 1 FOR iVal IN [0..(bnoise.dim-1)]: FOR jVal IN [0..(bnoise.dim-1)]: bnoise.map[iVal][jVal] = PRN() FOR iSurf in [0..1]: FOR iCtrl in [0..(nbCtrl-1)]: FOR jCtrl in [0..(nbCtrl-1)]: IF (iCtrl MODULO 3) != 0 AND (jCtrl MODULO 3) != 0: bnoise.surf[iSurf][iCtrl][jCtrl] = PRN() * order FOR iCtrl in [0..(nbCtrl-1)]: FOR jCtrl in [0..(nbCtrl-1)]: IF (iCtrl MODULO 3) != 0 AND (jCtrl MODULO 3) == 0: jCtrlPrev = (jCtrl - 1) MODULO (nbCtrl - 1) jCtrlNext = (jCtrl + 1) MODULO (nbCtrl - 1) bnoise.surf[iSurf][iCtrl][jCtrl] = 0.5 * (bnoise.surf[iSurf][iCtrl][jCtrlPrev] + bnoise.surf[iSurf][iCtrl][jCtrlNext]) ELSE IF (iCtrl MODULO 3) == 0 AND (jCtrl MODULO 3) != 0: iCtrlPrev = (iCtrl - 1) MODULO (nbCtrl - 1) iCtrlNext = (iCtrl + 1) MODULO (nbCtrl - 1) bnoise.surf[iSurf][iCtrl][jCtrl] = 0.5 * (bnoise.surf[iSurf][iCtrlPrev][jCtrl] + bnoise.surf[iSurf][iCtrlNext][jCtrl]) FOR iCtrl in [0..(nbCtrl-1)]: FOR jCtrl in [0..(nbCtrl-1)]: IF (iCtrl MODULO 3) == 0 AND (jCtrl MODULO 3) == 0: iCtrlPrev = (iCtrl - 1) MODULO (nbCtrl - 1) iCtrlNext = (iCtrl + 1) MODULO (nbCtrl - 1) bnoise.surf[iSurf][iCtrl][jCtrl] = 0.5 * (bnoise.surf[iSurf][iCtrlPrev][jCtrl] + bnoise.surf[iSurf][iCtrlNext][jCtrl]) FOR iCtrl in [0..(nbCtrl-1)]: FOR jCtrl in [0..2]: bnoise.surf[iSurf][iCtrl][nbCtrl + jCtrl] = bnoise.surf[iSurf][iCtrl][1 + jCtrl]; FOR iCtrl in [0..2]: FOR jCtrl in [0..(nbCtrl-1)]: bnoise.surf[iSurf][nbCtrl + iCtrl][jCtrl] = bnoise.surf[iSurf][1 + iCtrl][jCtrl]; FOR iCtrl in [0..2]: FOR jCtrl in [0..2]: bnoise.surf[iSurf][nbCtrl + iCtrl][nbCtrl + jCtrl] = bnoise.surf[iSurf][1 + iCtrl][1 + jCtrl]; FUNCTION Lerp( t: float, v: array of two floats): RETURN v[0] + t * (v[1] - v[0]) FUNCTION EvalBezier( t: float in [0,1], position on the curve, order: int, order of the curve curve: array of (order+1) float, the curve control): v: array of (order+1) floats v = copy of curve FOR depth IN [0..(order-1)]: FOR i IN [0..(order-depth-1)]: v[i] = Lerp(t, v[i..(i+1)]) RETURN v[0] FUNCTION EvalBezierSurf: t: array of 2 floats in [0,1], position on the surface, order: int, order of the surface surf: array of (order+1)*(order+1) floats, the surface controls): v: array of (order+1) floats FOR i IN [0..order]: v[i] = EvalBezier(t[0], surf[i], order) return EvalBezier(t[1], v, order); FUNCTION SmootherStep(x: float, the value to smooth): RETURN x * x * x * (x * (x * 6 - 15) + 10) FUNCTION BNoiseGet( bnoise: BNoise, the noise evaluated, pos: array of two floats in [-inf,inf], the input position, one tile is contained in [0,1]x[0,1]): modPos: array of two floats t: array of two floats idxCurve: array of two ints FOR iCoord IN [0..1]: modPos[iCoord] = bnoise.dim * (pos[iCoord] MODULO 1) t[iCoord] = modPos[iCoord] MODULO 1 idxCurve[iCoord] = FLOOR(modPos[iCoord]) mappedPos: array of two floats FOR iCoord IN [0..1]: surf = bnoise.surf[iCoord] [idxCurve[0]..(idxCurve[0]+3)][idxCurve[1]..(idxCurve[1]+3)] mappedPos[iCoord] = EvalBezierSurf(t, surf, 3) posIdx: array of 2 ints posCell: array of 2 floats FOR iCoord IN [0..1]: posIdx[iCoord] = FLOOR(mappedPos[iCoord]) posCell[iCoord] = SmootherStep(mappedPos[iCoord] MODULO 1) val: array of 2 floats FOR iCoord IN [0..1]: idx = (posIdx[1] + iCoord) * bnoise.dim + posIdx[0] val[iCoord] = Lerp(posCell[0], bnoise.map[idx..(idx+1)]) RETURN Lerp(posCell[1], val)

Edit on 2023/04/03: Note that the output value range for one instance of B-Noise generally won't cover the whole [0,1] range. This range is included in the grid values range, these values being picked at random they will most often not contain 0 and 1, hence their range, and consequently the one of the output, will be smaller than [0,1]. The lower the order of the noise, the higher the probablity for a small range of output values. If this is a concern for the user, normalising the grid values right after their initialisation will almost solve the problem. Only 'almost' because the values in the grid are interpolated at positions generated by the Bezier surfaces, and these positions are not guaranteed to cover the whole grid. That would probably be good enough for artistic purposes.

The examples below show the new results. We can clearly see that the directional bias has disappeared.

One tile, order 1, image #1 One tile, order 1, image #2 One tile, order 1, image #3 One tile, order 1, image #4 One tile, order 1, image #5 One tile, order 1, image #6

One tile, order 2, image #1 One tile, order 2, image #2 One tile, order 2, image #3 One tile, order 2, image #4 One tile, order 2, image #5 One tile, order 2, image #6

One tile, order 4, image #1 One tile, order 4, image #2 One tile, order 4, image #3 One tile, order 4, image #4 One tile, order 4, image #5 One tile, order 4, image #6

One tile, order 8, image #1 One tile, order 8, image #2 One tile, order 8, image #3 One tile, order 8, image #4 One tile, order 8, image #5 One tile, order 8, image #6

One tile, order 16, image #1 One tile, order 16, image #2 One tile, order 16, image #3 One tile, order 16, image #4 One tile, order 16, image #5 One tile, order 16, image #6

3x3 tiles, order 8, image #1 3x3 tiles, order 8, image #2 3x3 tiles, order 8, image #3 3x3 tiles, order 8, image #4 3x3 tiles, order 8, image #5 3x3 tiles, order 8, image #6

depth 5, lacunarity 2, gain 0.5, image #1 depth 5, lacunarity 2, gain 0.5, image #2 depth 5, lacunarity 2, gain 0.5, image #3 depth 5, lacunarity 2, gain 0.5, image #4

And finally, the implementation in C. Header file:

/* BNoise - A 2D noise generator. Copyright (C) 2023 Pascal Baillehache baillehache.pascal@gmail.com https://baillehachepascal.dev This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. */ #ifndef BNOISE_H #define BNOISE_H #include <stdlib.h> #include <string.h> #include <math.h> // BNoise structure typedef struct { // Dimension of the grid (square) in number of vertices int dim; float fDim; // Noise data float* map; float* surf[2]; int stride; } BNoise; // Allocate memory and create a BNoise instance // Inputs: // order: order of the noise. The grid will be of size (dim+1)*(dim+1) in // number of vertices and the Bezier curve will be of order (4+dim). // Output: // Return the new BNoise instance BNoise* BNoiseAlloc(int const order); // Free memory used by a BNoise // Inputs: // that: the BNoise void BNoiseFree(BNoise** const that); // Get the noise value at a given position // Inputs: // that: the BNoise // pos: the 2D position (one tile is contained in [0,1]x[0,1]) // Output: // Return the noise value at the position, in [0,1] float BNoiseGet( BNoise* const that, float const pos[2]); #endif

Body:

/* BNoise - A 2D noise generator. Copyright (C) 2023 Pascal Baillehache baillehache.pascal@gmail.com https://baillehachepascal.dev This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "bnoise.h" // Shortcut equivalent to a for loop iterating from 0 to (nbIter-1) // and storing the current iteration index in varName #define loop(varName, nbIter) \ for(__typeof__((nbIter) + 0) (varName) = 0; \ (varName) < (nbIter); ++(varName)) // Get a pseudo random number in [0,1] // Output: // Return the pseudo random number. static float rng(void) { return (float)rand() / (float)RAND_MAX; } // Allocate memory and create a BNoise instance // Inputs: // order: order of the noise. The grid will be of size (dim+1)*(dim+1) in // number of vertices and the Bezier curve will be of order (4+dim). // Output: // Return the new BNoise instance BNoise* BNoiseAlloc(int const order) { // Allocate memory BNoise* that = malloc(sizeof(BNoise)); // Set the dimension of the grid that->dim = order + 1; that->fDim = (float)(that->dim); // Allocate memory for the vertices value and initialise them to random // number in [0,1] size_t nbVertex = that->dim * that->dim; that->map = malloc(nbVertex * sizeof(float)); loop(i, nbVertex) that->map[i] = rng(); // Calculate the number of controls per axis in a surface int nbCtrl = that->dim * 3 + 1; // Calculate the stride, used to calculate indices of controls in the Bezier // surfaces, add 3 controls for the extra segment (cf below) that->stride = nbCtrl + 3; // Loop on the two Bezier surfaces loop(iSurf, 2) { // Allocate memory for the controls that->surf[iSurf] = malloc(sizeof(float) * that->stride * that->stride); // The center handles are set to random numbers in [0, order] float* ctrls = that->surf[iSurf]; loop(iCtrl, nbCtrl) loop(jCtrl, nbCtrl) { if((iCtrl % 3) != 0 && (jCtrl % 3) != 0) { ctrls[jCtrl * that->stride + iCtrl] = rng() * (float)order; } } // The other handles are set to the average of their neighbour handles loop(iCtrl, nbCtrl) loop(jCtrl, nbCtrl) { if((iCtrl % 3) != 0 && (jCtrl % 3) == 0) { int jCtrlPrev = (jCtrl - 1) % (nbCtrl - 1); if(jCtrlPrev < 0) jCtrlPrev += nbCtrl - 1; int jCtrlNext = (jCtrl + 1) % (nbCtrl - 1); ctrls[jCtrl * that->stride + iCtrl] = 0.5 * ( ctrls[jCtrlPrev * that->stride + iCtrl] + ctrls[jCtrlNext * that->stride + iCtrl] ); } else if((iCtrl % 3) == 0 && (jCtrl % 3) != 0) { int iCtrlPrev = (iCtrl - 1) % (nbCtrl - 1); if(iCtrlPrev < 0) iCtrlPrev += nbCtrl - 1; int iCtrlNext = (iCtrl + 1) % (nbCtrl - 1); ctrls[jCtrl * that->stride + iCtrl] = 0.5 * ( ctrls[jCtrl * that->stride + iCtrlPrev] + ctrls[jCtrl * that->stride + iCtrlNext] ); } } // Finally the anchors are set to the average of their neighbour handles loop(iCtrl, nbCtrl) loop(jCtrl, nbCtrl) { if((iCtrl % 3) == 0 && (jCtrl % 3) == 0) { int iCtrlPrev = (iCtrl - 1) % (nbCtrl - 1); if(iCtrlPrev < 0) iCtrlPrev += nbCtrl - 1; int iCtrlNext = (iCtrl + 1) % (nbCtrl - 1); ctrls[jCtrl * that->stride + iCtrl] = 0.5 * ( ctrls[jCtrl * that->stride + iCtrlPrev] + ctrls[jCtrl * that->stride + iCtrlNext] ); } } // Duplicate the first segments into extra segments to avoid a modulo in // BNoiseGet() loop(jCtrl, 3) loop(iCtrl, nbCtrl) { ctrls[(nbCtrl + jCtrl) * that->stride + iCtrl] = ctrls[(1 + jCtrl) * that->stride + iCtrl]; } loop(jCtrl, nbCtrl) loop(iCtrl, 3) { ctrls[jCtrl * that->stride + (iCtrl + nbCtrl)] = ctrls[jCtrl * that->stride + (iCtrl + 1)]; } loop(iCtrl, 3) loop(jCtrl, 3) { ctrls[(nbCtrl + jCtrl) * that->stride + (iCtrl + nbCtrl)] = ctrls[(1 + jCtrl) * that->stride + (iCtrl + 1)]; } } // Return the new noise instance return that; } // Free memory used by a BNoise // Inputs: // that: the BNoise void BNoiseFree(BNoise** const that) { if(that == NULL || *that == NULL) return; loop(iSurf, 2) free((*that)->surf[iSurf]); free((*that)->map); free(*that); *that = NULL; } // Smoother step function // Inputs: // x: the input value in [0,1] // Output: // Return the smoothed value, in [0, 1] static float SmootherStep(float const x) { return x * x * x * (x * (x * 6.0 - 15.0) + 10.0); } // LERP function, map a float value linearly from [0,1] to a range // Inputs: // x: the input // to: range of the output // Output: // Return the mapped input static float Lerp( float const x, float const* const to) { return to[0] + x * (to[1] - to[0]); } // Calculate the value of a Bezier curve // Inputs: // t: argument of the function, in [0, 1] // params: the control values of the Bezier // order: the order of the Bezier // Output: // Return the value of the Bezier static float EvalBezier( float const t, float const* const params, int const order) { float v[order + 1]; memcpy(v, params, sizeof(float) * (order + 1)); loop(depth, order) loop(i, order - depth) v[i] = Lerp(t, v + i); return v[0]; } // Calculate the value of a Bezier surface // Inputs: // t: argument of the function, in [0, 1] // params: the control values of the Bezier // stride: stride for the second dimension in control values // order: the order of the Bezier // Output: // Return the value of the Bezier static float EvalBezierSurf( float* const t, float const* const params, int const stride, int const order) { float v[order + 1]; loop(i, order + 1) v[i] = EvalBezier(t[0], params + i * stride, order); return EvalBezier(t[1], v, order); } // Get the noise value at a given position // Inputs: // that: the BNoise // pos: the 2D position (one tile is contained in [0,1]x[0,1]) // Output: // Return the noise value at the position, in [0,1] float BNoiseGet( BNoise* const that, float const pos[2]) { // Warp and scale the input position from [0, 1] over the interval // [0, that->dim] float modPos[2] = {0.0}; float t[2] = {0.0}; int idxCurve[2] = {0}; loop(iCoord, 2) { if(pos[iCoord] >= 0.0) { modPos[iCoord] = that->fDim * fmodf(pos[iCoord], 1.0); } else { modPos[iCoord] = that->fDim * (1.0 - fmodf(fabs(pos[iCoord]), 1.0)); } t[iCoord] = fmodf(modPos[iCoord], 1.0); idxCurve[iCoord] = floor(modPos[iCoord]); } // Map the warped and scaled position to a grid position using the two // Bezier surfaces float mappedPos[2] = {0.0}; loop(iCoord, 2) { mappedPos[iCoord] = EvalBezierSurf( t, that->surf[iCoord] + (idxCurve[1] * 3 * that->stride + idxCurve[0] * 3), that->stride, 3); } // Get the indices of the cell and local position in that cell for the // mapped position int posIdx[2] = {0}; float posCell[2] = {0.0}; loop(iCoord, 2) { posIdx[iCoord] = (int)floor(mappedPos[iCoord]); posCell[iCoord] = fmodf(mappedPos[iCoord], 1.0); // Apply smootherstep to the local position for derivability posCell[iCoord] = SmootherStep(posCell[iCoord]); } // Calculate the result value using bilinear interpolation on the cell // corners value using the local position float val[2] = {0.0}; loop(i, 2) { val[i] = Lerp(posCell[0], that->map + (posIdx[1] + i) * that->dim + posIdx[0]); } return Lerp(posCell[1], val); }

Download the source here. Any question or comment is welcome by email. The implementation in Javascript is available for download here, and the interactive demo has been updated to use the new version introduced in this article.

Edit on 2023/04/06: Correction of a bug in the C implementation. Initialisation of center handles used round(rng() * (float)order) instead of rng() * (float)order.