The advection vector doesn't have to be a constant: with some tweaking we can localize must of the vorticity to a single point on the sphre
Determining the distance to a point
The flow noise presented in a previous article has an advection parameter that determines how much small details are swept along with turbulent larger features. There is however nothing that prevents us to make this parameter dependent on the distance to some location to introduce some variation in the overall swirl patterns on the the sphere. An example is shown in the next video sample where the vorticity has its maximum near a point on the lower right front side:
Example node setup
The node setup looks like this:
The nodes in the yellow box calculate the distance to some hotspot ([1, -0.3, 1] not visible in this screenshot) and then calculate the lenght of this vector by determining the inproduct with itself. Then we add one to prevent division by zero in the next step, where we calculate the inverse. By raising it to some power we can control the 'tightness' of the spot, i.e. how fast our value decreases when we are further away from the hotspot. The final multiply node gives us some control on how strong the overall effect is.
In this article I present an OSL implementation of flow noise, both in 2D and in 4D
Flow Noise
in two previous posts (I, II) I showed a port to OSL of Stefan Gustavsons implementations of Simplex Noise and a fBM variant. Four dimensional Simplex Noise may be used to simulate the time evolution of 3D noise and in some situations like for example boiling liquids the result might be perfectly adequate. However if we want to capture the effects of swirls as seen in flowing liquids and gasses a better option is to use flow noise.
Flow noise was pioneered by Perlin and Neyret and the concept of adding vorticity (swirls) by rotating the lattice gradient vectors used to create the noise and adding higher octaves of noise at positions displaced by the lower octaves (advection) isn't limited to Perlin noise but can be applied to any lattice gradient noise, including Simplex Noise.
The code presented below is a straightforward and unoptimised port of Stefan Gustavsons original implementation in C. 2D and 3D versions of this noise are provided in a single shader,the type can be selected with the Dim parameter. By keyframing the Angle parameter you can create a swirling motion. The character of this motion may be influenced by altering the Octaves, lacunarity, H and advection parameters. The advection parameter basically determines how much small details in the flow will be swept along by the larger features.
Example node setup
The video of the swirling sphere was created with the following node setup:
The angle parameter was varied from 0 to 6.28 in 150 frames.
An additional trick to restrict the advection to some specific point is illustrated in another post.
The shader is an unwieldy large piece of code but nevertheless presented here in its entirety. Just click view plain to download it.
/* Simplex flow noise in 2 and 3D
*
* Based on Stefan Gustavsons orignal implementation in srdnoise23.c
* (http://webstaff.itn.liu.se/~stegu/simplexnoise/DSOnoises.html)
*
* OSL port Michel Anders (varkenvarken) 2013-02-04
* original comment is is left mostly in place, OSL specific comments
* are preceded by MJA
*
* This code was placed in the public domain by its original author,
* Stefan Gustavson. You may use it as you see fit, but
* attribution is appreciated.
*
*/
int FASTFLOOR(float x) {
int xi = (int)x;
return x < xi ? xi-1 : xi;
}
shader flownoise(
point Pos = P,
float Scale = 1,
int Octaves = 1,
float H = 1,
float lacunarity=2,
float Angle=0,
float advection=0,
int Dim = 2,
output float fac = 0
){
int perm[512] = {151,160,137,91,90,15,
131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
151,160,137,91,90,15,
131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180};
// MJA precomputing this table instead of calculating it as done in the
// original code saves 30% running time.
int permMod12[512] = { 7, 4, 5, 7, 6, 3, 11, 1, 9, 11, 0, 5, 2, 5, 7, 9, 8,
0, 7, 6, 9, 10, 8, 3, 1, 0, 9, 10, 11, 10, 6, 4, 7, 0, 6, 3, 0, 2, 5, 2, 10,
0, 3, 11, 9, 11, 11, 8, 9, 9, 9, 4, 9, 5, 8, 3, 6, 8, 5, 4, 3, 0, 8, 7, 2, 9,
11, 2, 7, 0, 3, 10, 5, 2, 2, 3, 11, 3, 1, 2, 0, 7, 1, 2, 4, 9, 8, 5, 7, 10,
5, 4, 4, 6, 11, 6, 5, 1, 3, 5, 1, 0, 8, 1, 5, 4, 0, 7, 4, 5, 6, 1, 8, 4, 3,
10, 8, 8, 3, 2, 8, 4, 1, 6, 5, 6, 3, 4, 4, 1, 10, 10, 4, 3, 5, 10, 2, 3, 10,
6, 3, 10, 1, 8, 3, 2, 11, 11, 11, 4, 10, 5, 2, 9, 4, 6, 7, 3, 2, 9, 11, 8, 8,
2, 8, 10, 7, 10, 5, 9, 5, 11, 11, 7, 4, 9, 9, 10, 3, 1, 7, 2, 0, 2, 7, 5, 8,
4, 10, 5, 4, 8, 2, 6, 1, 0, 11, 10, 2, 1, 10, 6, 0, 0, 11, 11, 6, 1, 9, 3, 1,
7, 9, 2, 11, 11, 1, 0, 10, 7, 1, 7, 10, 1, 4, 0, 0, 8, 7, 1, 2, 9, 7, 4, 6, 2,
6, 8, 1, 9, 6, 6, 7, 5, 0, 0, 3, 9, 8, 3, 6, 6, 11, 1, 0, 0, 7, 4, 5, 7, 6, 3,
11, 1, 9, 11, 0, 5, 2, 5, 7, 9, 8, 0, 7, 6, 9, 10, 8, 3, 1, 0, 9, 10, 11, 10,
6, 4, 7, 0, 6, 3, 0, 2, 5, 2, 10, 0, 3, 11, 9, 11, 11, 8, 9, 9, 9, 4, 9, 5, 8,
3, 6, 8, 5, 4, 3, 0, 8, 7, 2, 9, 11, 2, 7, 0, 3, 10, 5, 2, 2, 3, 11, 3, 1, 2,
0, 7, 1, 2, 4, 9, 8, 5, 7, 10, 5, 4, 4, 6, 11, 6, 5, 1, 3, 5, 1, 0, 8, 1, 5, 4,
0, 7, 4, 5, 6, 1, 8, 4, 3, 10, 8, 8, 3, 2, 8, 4, 1, 6, 5, 6, 3, 4, 4, 1, 10, 10,
4, 3, 5, 10, 2, 3, 10, 6, 3, 10, 1, 8, 3, 2, 11, 11, 11, 4, 10, 5, 2, 9, 4, 6, 7,
3, 2, 9, 11, 8, 8, 2, 8, 10, 7, 10, 5, 9, 5, 11, 11, 7, 4, 9, 9, 10, 3, 1, 7, 2,
0, 2, 7, 5, 8, 4, 10, 5, 4, 8, 2, 6, 1, 0, 11, 10, 2, 1, 10, 6, 0, 0, 11, 11, 6,
1, 9, 3, 1, 7, 9, 2, 11, 11, 1, 0, 10, 7, 1, 7, 10, 1, 4, 0, 0, 8, 7, 1, 2, 9, 7,
4, 6, 2, 6, 8, 1, 9, 6, 6, 7, 5, 0, 0, 3, 9, 8, 3, 6, 6, 11, 1, 0, 0};
/*
* Gradient tables. These could be programmed the Ken Perlin way with
* some clever bit-twiddling, but this is more clear, and not really slower.
*/
vector grad2[8] = {
vector( -1.0, -1.0 , 0.0), vector( 1.0, 0.0 , 0.0) , vector( -1.0, 0.0 , 0.0) , vector( 1.0, 1.0 , 0.0) ,
vector( -1.0, 1.0 , 0.0) , vector( 0.0, -1.0 , 0.0) , vector( 0.0, 1.0 , 0.0) , vector( 1.0, -1.0 , 0.0)
};
/*
* For 3D, we define two orthogonal vectors in the desired rotation plane.
* These vectors are based on the midpoints of the 12 edges of a cube,
* they all rotate in their own plane and are never coincident or collinear.
* A larger array of random vectors would also do the job, but these 12
* (including 4 repeats to make the array length a power of two) work better.
* They are not random, they are carefully chosen to represent a small
* isotropic set of directions for any rotation angle.
*/
/* a = sqrt(2)/sqrt(3) = 0.816496580 */
#define a 0.81649658
vector grad3u[16] = {
vector( 1.0, 0.0, 1.0 ), vector( 0.0, 1.0, 1.0 ), // 12 cube edges
vector( -1.0, 0.0, 1.0 ), vector( 0.0, -1.0, 1.0 ),
vector( 1.0, 0.0, -1.0 ), vector( 0.0, 1.0, -1.0 ),
vector( -1.0, 0.0, -1.0 ), vector( 0.0, -1.0, -1.0 ),
vector( a, a, a ), vector( -a, a, -a ),
vector( -a, -a, a ), vector( a, -a, -a ),
vector( -a, a, a ), vector( a, -a, a ),
vector( a, -a, -a ), vector( -a, a, -a )
};
vector grad3v[16] = {
vector( -a, a, a ), vector( -a, -a, a ),
vector( a, -a, a ), vector( a, a, a ),
vector( -a, -a, -a ), vector( a, -a, -a ),
vector( a, a, -a ), vector( -a, a, -a ),
vector( 1.0, -1.0, 0.0 ), vector( 1.0, 1.0, 0.0 ),
vector( -1.0, 1.0, 0.0 ), vector( -1.0, -1.0, 0.0 ),
vector( 1.0, 0.0, 1.0 ), vector( -1.0, 0.0, 1.0 ), // 4 repeats to make 16
vector( 0.0, 1.0, -1.0 ), vector( 0.0, -1.0, -1.0 )
};
#undef a
/*
* Helper functions to compute rotated gradients and
* gradients-dot-residualvectors in 2D and 3D.
*/
void gradrot2( int hash, float sin_t, float cos_t, float gx, float gy ) {
int h = hash & 7;
float gx0 = grad2[h][0];
float gy0 = grad2[h][1];
gx = cos_t * gx0 - sin_t * gy0;
gy = sin_t * gx0 + cos_t * gy0;
return;
}
void gradrot3( int hash, float sin_t, float cos_t, float gx, float gy, float gz ) {
int h = hash & 15;
float gux = grad3u[h][0];
float guy = grad3u[h][1];
float guz = grad3u[h][2];
float gvx = grad3v[h][0];
float gvy = grad3v[h][1];
float gvz = grad3v[h][2];
gx = cos_t * gux + sin_t * gvx;
gy = cos_t * guy + sin_t * gvy;
gz = cos_t * guz + sin_t * gvz;
return;
}
float graddotp3( float gx, float gy, float gz, float x, float y, float z ) {
return gx * x + gy * y + gz * z;
}
float angle = Angle; // MJA copy input parameter so that we can write it
float pwr = 1.0;
float pwHL = pow(lacunarity,-H);
/* Skewing factors for 2D simplex grid:
* F2 = 0.5*(sqrt(3.0)-1.0)
* G2 = (3.0-Math.sqrt(3.0))/6.0
*/
#define F2 0.366025403
#define G2 0.211324865
/** 2D simplex noise with derivatives.
* If the last two arguments are not null, the analytic derivative
* (the 2D gradient of the scalar noise field) is also calculated.
*/
if( Dim == 2 ){
float x = Pos[0]*Scale;
float y = Pos[1]*Scale;
float dnoise_dx;
float dnoise_dy;
float sin_t, cos_t; /* Sine and cosine for the gradient rotation angle */
sin_t = sin( angle );
cos_t = cos( angle );
for(int p=0; p < Octaves; p++){
float n0, n1, n2; /* Noise contributions from the three simplex corners */
float gx0, gy0, gx1, gy1, gx2, gy2; /* Gradients at simplex corners */
/* Skew the input space to determine which simplex cell we're in */
float s = ( x + y ) * F2; /* Hairy factor for 2D */
float xs = x + s;
float ys = y + s;
int i = FASTFLOOR( xs );
int j = FASTFLOOR( ys );
float t = ( i + j ) * G2;
float X0 = i - t; /* Unskew the cell origin back to (x,y) space */
float Y0 = j - t;
float x0 = x - X0; /* The x,y distances from the cell origin */
float y0 = y - Y0;
/* For the 2D case, the simplex shape is an equilateral triangle.
* Determine which simplex we are in. */
int i1, j1; /* Offsets for second (middle) corner of simplex in (i,j) coords */
if( x0 > y0 ) { i1 = 1; j1 = 0; } /* lower triangle, XY order: (0,0)->(1,0)->(1,1) */
else { i1 = 0; j1 = 1; } /* upper triangle, YX order: (0,0)->(0,1)->(1,1) */
/* A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
* a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
* c = (3-sqrt(3))/6 */
float x1 = x0 - i1 + G2; /* Offsets for middle corner in (x,y) unskewed coords */
float y1 = y0 - j1 + G2;
float x2 = x0 - 1.0 + 2.0 * G2; /* Offsets for last corner in (x,y) unskewed coords */
float y2 = y0 - 1.0 + 2.0 * G2;
/* Wrap the integer indices at 256, to avoid indexing perm[] out of bounds */
int ii = i & 255; // MJA was % 256 but OSL mod is not the same as C %
int jj = j & 255;
/* Calculate the contribution from the three corners */
float t0 = 0.5 - x0 * x0 - y0 * y0;
float t20, t40;
if( t0 < 0.0 ) t40 = t20 = t0 = n0 = gx0 = gy0 = 0.0; /* No influence */
else {
gradrot2( perm[ii + perm[jj]], sin_t, cos_t, gx0, gy0 );
t20 = t0 * t0;
t40 = t20 * t20;
n0 = t40 * ( gx0 * x0 + gy0 * y0 );
}
float t1 = 0.5 - x1 * x1 - y1 * y1;
float t21, t41;
if( t1 < 0.0 ) t21 = t41 = t1 = n1 = gx1 = gy1 = 0.0; /* No influence */
else {
gradrot2( perm[ii + i1 + perm[jj + j1]], sin_t, cos_t, gx1, gy1 );
t21 = t1 * t1;
t41 = t21 * t21;
n1 = t41 * ( gx1 * x1 + gy1 * y1 );
}
float t2 = 0.5 - x2 * x2 - y2 * y2;
float t22, t42;
if( t2 < 0.0 ) t42 = t22 = t2 = n2 = gx2 = gy2 = 0.0; /* No influence */
else {
gradrot2( perm[ii + 1 + perm[jj + 1]], sin_t, cos_t, gx2, gy2 );
t22 = t2 * t2;
t42 = t22 * t22;
n2 = t42 * ( gx2 * x2 + gy2 * y2 );
}
/* Add contributions from each corner to get the final noise value.
* The result is scaled to return values in the interval [-1,1]. */
float noise = 70.0 * ( n0 + n1 + n2 ); // MJA scale factor was 40
/* Compute derivative, if requested by supplying non-null pointers
* for the last two arguments */
if( advection != 0 ){
/* A straight, unoptimised calculation would be like:
* *dnoise_dx = -8.0 * t20 * t0 * x0 * ( gx0 * x0 + gy0 * y0 ) + t40 * gx0;
* *dnoise_dy = -8.0 * t20 * t0 * y0 * ( gx0 * x0 + gy0 * y0 ) + t40 * gy0;
* *dnoise_dx += -8.0 * t21 * t1 * x1 * ( gx1 * x1 + gy1 * y1 ) + t41 * gx1;
* *dnoise_dy += -8.0 * t21 * t1 * y1 * ( gx1 * x1 + gy1 * y1 ) + t41 * gy1;
* *dnoise_dx += -8.0 * t22 * t2 * x2 * ( gx2 * x2 + gy2 * y2 ) + t42 * gx2;
* *dnoise_dy += -8.0 * t22 * t2 * y2 * ( gx2 * x2 + gy2 * y2 ) + t42 * gy2;
*/
float temp0 = t20 * t0 * ( gx0* x0 + gy0 * y0 );
dnoise_dx = temp0 * x0;
dnoise_dy = temp0 * y0;
float temp1 = t21 * t1 * ( gx1 * x1 + gy1 * y1 );
dnoise_dx += temp1 * x1;
dnoise_dy += temp1 * y1;
float temp2 = t22 * t2 * ( gx2* x2 + gy2 * y2 );
dnoise_dx += temp2 * x2;
dnoise_dy += temp2 * y2;
dnoise_dx *= -8.0;
dnoise_dy *= -8.0;
dnoise_dx += t40 * gx0 + t41 * gx1 + t42 * gx2;
dnoise_dy += t40 * gy0 + t41 * gy1 + t42 * gy2;
dnoise_dx *= 70.0; /* Scale derivative to match the noise scaling */ // MJA scale factor was 40
dnoise_dy *= 70.0;
float advect_x = -advection * dnoise_dx;
float advect_y = -advection * dnoise_dy;
fac += noise * pwr;
pwr *= pwHL;
x *= lacunarity;
y *= lacunarity;
x += advect_x*pow(lacunarity,p);
y += advect_y*pow(lacunarity,p);
angle *= lacunarity;
}else{
fac += noise * pwr;
pwr *= pwHL;
x *= lacunarity;
y *= lacunarity;
}
}
}else if(Dim==3){
/* Skewing factors for 3D simplex grid:
* F3 = 1/3
* G3 = 1/6 */
#define F3 0.333333333
#define G3 0.166666667
float x=Pos[0]*Scale;
float y=Pos[1]*Scale;
float z=Pos[2]*Scale;
float dnoise_dx;
float dnoise_dy;
float dnoise_dz;
float n0, n1, n2, n3; /* Noise contributions from the four simplex corners */
float noise; /* Return value */
float gx0, gy0, gz0, gx1, gy1, gz1; /* Gradients at simplex corners */
float gx2, gy2, gz2, gx3, gy3, gz3;
float sin_t, cos_t; /* Sine and cosine for the gradient rotation angle */
sin_t = sin( angle );
cos_t = cos( angle );
for(int p=0; p < Octaves; p++){
/* Skew the input space to determine which simplex cell we're in */
float s = (x+y+z)*F3; /* Very nice and simple skew factor for 3D */
float xs = x+s;
float ys = y+s;
float zs = z+s;
int i = FASTFLOOR(xs);
int j = FASTFLOOR(ys);
int k = FASTFLOOR(zs);
float t = (float)(i+j+k)*G3;
float X0 = i-t; /* Unskew the cell origin back to (x,y,z) space */
float Y0 = j-t;
float Z0 = k-t;
float x0 = x-X0; /* The x,y,z distances from the cell origin */
float y0 = y-Y0;
float z0 = z-Z0;
/* For the 3D case, the simplex shape is a slightly irregular tetrahedron.
* Determine which simplex we are in. */
int i1, j1, k1; /* Offsets for second corner of simplex in (i,j,k) coords */
int i2, j2, k2; /* Offsets for third corner of simplex in (i,j,k) coords */
/* TODO: This code would benefit from a backport from the GLSL version! */
if(x0>=y0) {
if(y0>=z0)
{ i1=1; j1=0; k1=0; i2=1; j2=1; k2=0; } /* X Y Z order */
else if(x0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=0; k2=1; } /* X Z Y order */
else { i1=0; j1=0; k1=1; i2=1; j2=0; k2=1; } /* Z X Y order */
}
else { // x0 < y0
if(y0 < z0) { i1=0; j1=0; k1=1; i2=0; j2=1; k2=1; } /* Z Y X order */
else if(x0 < z0) { i1=0; j1=1; k1=0; i2=0; j2=1; k2=1; } /* Y Z X order */
else { i1=0; j1=1; k1=0; i2=1; j2=1; k2=0; } /* Y X Z order */
}
/* A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z),
* a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and
* a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where
* c = 1/6. */
float x1 = x0 - i1 + G3; /* Offsets for second corner in (x,y,z) coords */
float y1 = y0 - j1 + G3;
float z1 = z0 - k1 + G3;
float x2 = x0 - i2 + 2.0 * G3; /* Offsets for third corner in (x,y,z) coords */
float y2 = y0 - j2 + 2.0 * G3;
float z2 = z0 - k2 + 2.0 * G3;
float x3 = x0 - 1.0 + 3.0 * G3; /* Offsets for last corner in (x,y,z) coords */
float y3 = y0 - 1.0 + 3.0 * G3;
float z3 = z0 - 1.0 + 3.0 * G3;
/* Wrap the integer indices at 256, to avoid indexing perm[] out of bounds */
int ii = i & 255; // MJA was % 256 but OSL mod is not the same as C %
int jj = j & 255;
int kk = k & 255;
/* Calculate the contribution from the four corners */
float t0 = 0.6 - x0*x0 - y0*y0 - z0*z0;
float t20, t40;
if(t0 < 0.0) n0 = t0 = t20 = t40 = gx0 = gy0 = gz0 = 0.0;
else {
gradrot3( perm[ii + perm[jj + perm[kk]]], sin_t, cos_t, gx0, gy0, gz0 );
t20 = t0 * t0;
t40 = t20 * t20;
n0 = t40 * graddotp3( gx0, gy0, gz0, x0, y0, z0 );
}
float t1 = 0.6 - x1*x1 - y1*y1 - z1*z1;
float t21, t41;
if(t1 < 0.0) n1 = t1 = t21 = t41 = gx1 = gy1 = gz1 = 0.0;
else {
gradrot3( perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]], sin_t, cos_t, gx1, gy1, gz1 );
t21 = t1 * t1;
t41 = t21 * t21;
n1 = t41 * graddotp3( gx1, gy1, gz1, x1, y1, z1 );
}
float t2 = 0.6 - x2*x2 - y2*y2 - z2*z2;
float t22, t42;
if(t2 < 0.0) n2 = t2 = t22 = t42 = gx2 = gy2 = gz2 = 0.0;
else {
gradrot3( perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]], sin_t, cos_t, gx2, gy2, gz2 );
t22 = t2 * t2;
t42 = t22 * t22;
n2 = t42 * graddotp3( gx2, gy2, gz2, x2, y2, z2 );
}
float t3 = 0.6 - x3*x3 - y3*y3 - z3*z3;
float t23, t43;
if(t3 < 0.0) n3 = t3 = t23 = t43 = gx3 = gy3 = gz3 = 0.0;
else {
gradrot3( perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]], sin_t, cos_t, gx3, gy3, gz3 );
t23 = t3 * t3;
t43 = t23 * t23;
n3 = t43 * graddotp3( gx3, gy3, gz3, x3, y3, z3 );
}
/* Add contributions from each corner to get the final noise value.
* The result is scaled to return values in the range [-1,1] */
noise = 28.0 * (n0 + n1 + n2 + n3);
/* Compute derivative, if requested by supplying non-null pointers
* for the last three arguments */
if(advection != 0){
/* A straight, unoptimised calculation would be like:
* *dnoise_dx = -8.0f * t20 * t0 * x0 * graddotp3(gx0, gy0, gz0, x0, y0, z0) + t40 * gx0;
* *dnoise_dy = -8.0f * t20 * t0 * y0 * graddotp3(gx0, gy0, gz0, x0, y0, z0) + t40 * gy0;
* *dnoise_dz = -8.0f * t20 * t0 * z0 * graddotp3(gx0, gy0, gz0, x0, y0, z0) + t40 * gz0;
* *dnoise_dx += -8.0f * t21 * t1 * x1 * graddotp3(gx1, gy1, gz1, x1, y1, z1) + t41 * gx1;
* *dnoise_dy += -8.0f * t21 * t1 * y1 * graddotp3(gx1, gy1, gz1, x1, y1, z1) + t41 * gy1;
* *dnoise_dz += -8.0f * t21 * t1 * z1 * graddotp3(gx1, gy1, gz1, x1, y1, z1) + t41 * gz1;
* *dnoise_dx += -8.0f * t22 * t2 * x2 * graddotp3(gx2, gy2, gz2, x2, y2, z2) + t42 * gx2;
* *dnoise_dy += -8.0f * t22 * t2 * y2 * graddotp3(gx2, gy2, gz2, x2, y2, z2) + t42 * gy2;
* *dnoise_dz += -8.0f * t22 * t2 * z2 * graddotp3(gx2, gy2, gz2, x2, y2, z2) + t42 * gz2;
* *dnoise_dx += -8.0f * t23 * t3 * x3 * graddotp3(gx3, gy3, gz3, x3, y3, z3) + t43 * gx3;
* *dnoise_dy += -8.0f * t23 * t3 * y3 * graddotp3(gx3, gy3, gz3, x3, y3, z3) + t43 * gy3;
* *dnoise_dz += -8.0f * t23 * t3 * z3 * graddotp3(gx3, gy3, gz3, x3, y3, z3) + t43 * gz3;
*/
float temp0 = t20 * t0 * graddotp3( gx0, gy0, gz0, x0, y0, z0 );
dnoise_dx = temp0 * x0;
dnoise_dy = temp0 * y0;
dnoise_dz = temp0 * z0;
float temp1 = t21 * t1 * graddotp3( gx1, gy1, gz1, x1, y1, z1 );
dnoise_dx += temp1 * x1;
dnoise_dy += temp1 * y1;
dnoise_dz += temp1 * z1;
float temp2 = t22 * t2 * graddotp3( gx2, gy2, gz2, x2, y2, z2 );
dnoise_dx += temp2 * x2;
dnoise_dy += temp2 * y2;
dnoise_dz += temp2 * z2;
float temp3 = t23 * t3 * graddotp3( gx3, gy3, gz3, x3, y3, z3 );
dnoise_dx += temp3 * x3;
dnoise_dy += temp3 * y3;
dnoise_dz += temp3 * z3;
dnoise_dx *= -8.0;
dnoise_dy *= -8.0;
dnoise_dz *= -8.0;
/* This corrects a bug in the original implementation */
dnoise_dx += t40 * gx0 + t41 * gx1 + t42 * gx2 + t43 * gx3;
dnoise_dy += t40 * gy0 + t41 * gy1 + t42 * gy2 + t43 * gy3;
dnoise_dz += t40 * gz0 + t41 * gz1 + t42 * gz2 + t43 * gz3;
dnoise_dx *= 28.0; /* Scale derivative to match the noise scaling */
dnoise_dy *= 28.0;
dnoise_dz *= 28.0;
float advect_x = -advection * dnoise_dx;
float advect_y = -advection * dnoise_dy;
float advect_z = -advection * dnoise_dz;
fac += noise * pwr;
pwr *= pwHL;
x *= lacunarity;
y *= lacunarity;
x += advect_x*pow(lacunarity,p);
y += advect_y*pow(lacunarity,p);
y += advect_z*pow(lacunarity,p);
angle *= lacunarity;
}else{
fac += noise * pwr;
pwr *= pwHL;
x *= lacunarity;
y *= lacunarity;
z *= lacunarity;
}
}
}
}
In this article a small update on my previous implementation of Simplex Noise for OSL: an extra feature to add turbulence in the form of fractional brownian motion (fBM) to the noise
Fractional Brownian Motion
In a previous article I showed how to implement 2, 3 and 4 dimensional simplex noise in OSL. In this article I show how this shader can be augmented with fractional brownian motion. Examples for 1 (= unperturbed noise), 2 and 3 octaves are shown here:
The code shown below isn't much altered from the previous version except for a loop around the core noise generation to add and scale additional octaves of noise. The shader node it creates gets extra input sockets for octaves, lacunarity and H which have the same meaning as in the regular Perlin implementation of fBM enhanced noise.
/*
* A speed-improved simplex noise algorithm for 2D, 3D and 4D in OSL.
*
* Based on example Java code by Stefan Gustavson (stegu@itn.liu.se).
* Optimisations by Peter Eastman (peastman@drizzle.stanford.edu).
* Better rank ordering method by Stefan Gustavson in 2012.
*
* This could be speeded up even further, but it's useful as it is.
*
* OSL port Michel Anders (varkenvarken) 2013-02-04
* original comment is is left in place, OSL specific comments
* are preceded by MJA
*
* This code was placed in the public domain by its original author,
* Stefan Gustavson. You may use it as you see fit, but
* attribution is appreciated.
*
* 2013-02-09 (varkenvarken) added fBM turbulence (Octaves)
*/
int fastfloor(float x) {
int xi = (int)x;
return x < xi ? xi-1 : xi;
}
// MJA it is safe to overload functions in OSL
// so here's some extra dot() versions
float dot(vector g, float x, float y) {
return g[0]*x + g[1]*y; }
float dot(vector g, float x, float y, float z) {
return g[0]*x + g[1]*y + g[2]*z; }
float dot(vector g, float t, float x, float y, float z, float w) {
return g[0]*x + g[1]*y + g[2]*z + t*w; }
shader simplexnoise(
point Pos = P,
float t = 0,
float Scale = 1,
int Octaves = 1,
float H = 1,
float lacunarity=2,
int Dim = 2,
output float fac = 0
){
vector grad3[12] = {vector(1,1,0),vector(-1,1,0),
vector(1,-1,0),vector(-1,-1,0),
vector(1,0,1),vector(-1,0,1),
vector(1,0,-1),vector(-1,0,-1),
vector(0,1,1),vector(0,-1,1),
vector(0,1,-1),vector(0,-1,-1)};
// MJA I couldn't get OSL to compile an array of structs so
// I separated a vec4 into a vector and a float
vector grad4v[32]= {vector(0,1,1) ,vector(0,1,1) ,
vector(0,1,-1) ,vector(0,1,-1),
vector(0,-1,1),vector(0,-1,1) ,
vector(0,-1,-1),vector(0,-1,-1),
vector(1,0,1) ,vector(1,0,1) ,
vector(1,0,-1) ,vector(1,0,-1),
vector(-1,0,1),vector(-1,0,1) ,
vector(-1,0,-1),vector(-1,0,-1),
vector(1,1,0) ,vector(1,1,0) ,
vector(1,-1,0) ,vector(1,-1,0),
vector(-1,1,0),vector(-1,1,0) ,
vector(-1,-1,0),vector(-1,-1,0),
vector(1,1,1) ,vector(1,1,-1) ,
vector(1,-1,1) ,vector(1,-1,-1),
vector(-1,1,1),vector(-1,1,-1),
vector(-1,-1,1),vector(-1,-1,-1)};
float grad4t[32]= { 1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,
1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,
0,0,0,0,0,0,0,0};
int perm[512] = {151,160,137,91,90,15,
131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
151,160,137,91,90,15,
131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180};
// MJA precomputing this table instead of calculating it as done in the
// original code saves 30% running time.
int permMod12[512] = { 7, 4, 5, 7, 6, 3, 11, 1, 9, 11, 0, 5, 2, 5, 7, 9, 8,
0, 7, 6, 9, 10, 8, 3, 1, 0, 9, 10, 11, 10, 6, 4, 7, 0, 6, 3, 0, 2, 5, 2, 10,
0, 3, 11, 9, 11, 11, 8, 9, 9, 9, 4, 9, 5, 8, 3, 6, 8, 5, 4, 3, 0, 8, 7, 2, 9,
11, 2, 7, 0, 3, 10, 5, 2, 2, 3, 11, 3, 1, 2, 0, 7, 1, 2, 4, 9, 8, 5, 7, 10,
5, 4, 4, 6, 11, 6, 5, 1, 3, 5, 1, 0, 8, 1, 5, 4, 0, 7, 4, 5, 6, 1, 8, 4, 3,
10, 8, 8, 3, 2, 8, 4, 1, 6, 5, 6, 3, 4, 4, 1, 10, 10, 4, 3, 5, 10, 2, 3, 10,
6, 3, 10, 1, 8, 3, 2, 11, 11, 11, 4, 10, 5, 2, 9, 4, 6, 7, 3, 2, 9, 11, 8, 8,
2, 8, 10, 7, 10, 5, 9, 5, 11, 11, 7, 4, 9, 9, 10, 3, 1, 7, 2, 0, 2, 7, 5, 8,
4, 10, 5, 4, 8, 2, 6, 1, 0, 11, 10, 2, 1, 10, 6, 0, 0, 11, 11, 6, 1, 9, 3, 1,
7, 9, 2, 11, 11, 1, 0, 10, 7, 1, 7, 10, 1, 4, 0, 0, 8, 7, 1, 2, 9, 7, 4, 6, 2,
6, 8, 1, 9, 6, 6, 7, 5, 0, 0, 3, 9, 8, 3, 6, 6, 11, 1, 0, 0, 7, 4, 5, 7, 6, 3,
11, 1, 9, 11, 0, 5, 2, 5, 7, 9, 8, 0, 7, 6, 9, 10, 8, 3, 1, 0, 9, 10, 11, 10,
6, 4, 7, 0, 6, 3, 0, 2, 5, 2, 10, 0, 3, 11, 9, 11, 11, 8, 9, 9, 9, 4, 9, 5, 8,
3, 6, 8, 5, 4, 3, 0, 8, 7, 2, 9, 11, 2, 7, 0, 3, 10, 5, 2, 2, 3, 11, 3, 1, 2,
0, 7, 1, 2, 4, 9, 8, 5, 7, 10, 5, 4, 4, 6, 11, 6, 5, 1, 3, 5, 1, 0, 8, 1, 5, 4,
0, 7, 4, 5, 6, 1, 8, 4, 3, 10, 8, 8, 3, 2, 8, 4, 1, 6, 5, 6, 3, 4, 4, 1, 10, 10,
4, 3, 5, 10, 2, 3, 10, 6, 3, 10, 1, 8, 3, 2, 11, 11, 11, 4, 10, 5, 2, 9, 4, 6, 7,
3, 2, 9, 11, 8, 8, 2, 8, 10, 7, 10, 5, 9, 5, 11, 11, 7, 4, 9, 9, 10, 3, 1, 7, 2,
0, 2, 7, 5, 8, 4, 10, 5, 4, 8, 2, 6, 1, 0, 11, 10, 2, 1, 10, 6, 0, 0, 11, 11, 6,
1, 9, 3, 1, 7, 9, 2, 11, 11, 1, 0, 10, 7, 1, 7, 10, 1, 4, 0, 0, 8, 7, 1, 2, 9, 7,
4, 6, 2, 6, 8, 1, 9, 6, 6, 7, 5, 0, 0, 3, 9, 8, 3, 6, 6, 11, 1, 0, 0};
// Skewing and unskewing factors for 2, 3, and 4 dimensions
float F2 = 0.5*(sqrt(3.0)-1.0);
float G2 = (3.0-sqrt(3.0))/6.0;
float F3 = 1.0/3.0;
float G3 = 1.0/6.0;
float F4 = (sqrt(5.0)-1.0)/4.0;
float G4 = (5.0-sqrt(5.0))/20.0;
float pwr = 1.0;
float pwHL = pow(lacunarity,-H);
if(Dim == 2){
// 2D simplex noise
fac = 0;
float xin=Pos[0]*Scale, yin=Pos[1]*Scale;
for(int p=0; p float n0, n1, n2; // Noise contributions from the three corners
// Skew the input space to determine which simplex cell we're in
float s = (xin+yin)*F2; // Hairy factor for 2D
int i = fastfloor(xin+s);
int j = fastfloor(yin+s);
float t = (i+j)*G2;
float X0 = i-t; // Unskew the cell origin back to (x,y) space
float Y0 = j-t;
float x0 = xin-X0; // The x,y distances from the cell origin
float y0 = yin-Y0;
// For the 2D case, the simplex shape is an equilateral triangle.
// Determine which simplex we are in.
int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords
if(x0>y0) {i1=1; j1=0;} // lower triangle, XY order: (0,0)->(1,0)->(1,1)
else {i1=0; j1=1;} // upper triangle, YX order: (0,0)->(0,1)->(1,1)
// A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
// a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
// c = (3-sqrt(3))/6
float x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords
float y1 = y0 - j1 + G2;
float x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords
float y2 = y0 - 1.0 + 2.0 * G2;
// Work out the hashed gradient indices of the three simplex corners
int ii = i & 255;
int jj = j & 255;
int gi0 = permMod12[ii+perm[jj]];
int gi1 = permMod12[ii+i1+perm[jj+j1]];
int gi2 = permMod12[ii+1+perm[jj+1]];
// Calculate the contribution from the three corners
float t0 = 0.5 - x0*x0-y0*y0;
if(t0 < 0) n0 = 0.0;
else {
t0 *= t0;
n0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient
}
float t1 = 0.5 - x1*x1-y1*y1;
if(t1 < 0) n1 = 0.0;
else {
t1 *= t1;
n1 = t1 * t1 * dot(grad3[gi1], x1, y1);
}
float t2 = 0.5 - x2*x2-y2*y2;
if(t2 < 0) n2 = 0.0;
else {
t2 *= t2;
n2 = t2 * t2 * dot(grad3[gi2], x2, y2);
}
// Add contributions from each corner to get the final noise value.
// The result is scaled to return values in the interval [-1,1].
fac += (70.0 * (n0 + n1 + n2)) * pwr;
pwr *= pwHL;
xin *= lacunarity;
yin *= lacunarity;
}
} else if(Dim == 3){
// 3D simplex noise
fac = 0 ;
float xin=Pos[0]*Scale, yin=Pos[1]*Scale, zin=Pos[2]*Scale;
for(int p=0; p float n0, n1, n2, n3; // Noise contributions from the four corners
// Skew the input space to determine which simplex cell we're in
float s = (xin+yin+zin)*F3; // Very nice and simple skew factor for 3D
int i = fastfloor(xin+s);
int j = fastfloor(yin+s);
int k = fastfloor(zin+s);
float t = (i+j+k)*G3;
float X0 = i-t; // Unskew the cell origin back to (x,y,z) space
float Y0 = j-t;
float Z0 = k-t;
float x0 = xin-X0; // The x,y,z distances from the cell origin
float y0 = yin-Y0;
float z0 = zin-Z0;
// For the 3D case, the simplex shape is a slightly irregular tetrahedron.
// Determine which simplex we are in.
int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords
int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords
if(x0>=y0) {
if(y0>=z0)
{ i1=1; j1=0; k1=0; i2=1; j2=1; k2=0; } // X Y Z order
else if(x0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=0; k2=1; } // X Z Y order
else { i1=0; j1=0; k1=1; i2=1; j2=0; k2=1; } // Z X Y order
}
else { // x0 < y0
if(y0 < z0) { i1=0; j1=0; k1=1; i2=0; j2=1; k2=1; } // Z Y X order
else if(x0 < z0) { i1=0; j1=1; k1=0; i2=0; j2=1; k2=1; } // Y Z X order
else { i1=0; j1=1; k1=0; i2=1; j2=1; k2=0; } // Y X Z order
}
// A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z),
// a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and
// a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where
// c = 1/6.
float x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z) coords
float y1 = y0 - j1 + G3;
float z1 = z0 - k1 + G3;
float x2 = x0 - i2 + 2.0*G3; // Offsets for third corner in (x,y,z) coords
float y2 = y0 - j2 + 2.0*G3;
float z2 = z0 - k2 + 2.0*G3;
float x3 = x0 - 1.0 + 3.0*G3; // Offsets for last corner in (x,y,z) coords
float y3 = y0 - 1.0 + 3.0*G3;
float z3 = z0 - 1.0 + 3.0*G3;
// Work out the hashed gradient indices of the four simplex corners
int ii = i & 255;
int jj = j & 255;
int kk = k & 255;
int gi0 = permMod12[ii+perm[jj+perm[kk]]];
int gi1 = permMod12[ii+i1+perm[jj+j1+perm[kk+k1]]];
int gi2 = permMod12[ii+i2+perm[jj+j2+perm[kk+k2]]];
int gi3 = permMod12[ii+1+perm[jj+1+perm[kk+1]]];
// Calculate the contribution from the four corners
float t0 = 0.6 - x0*x0 - y0*y0 - z0*z0;
if(t0 < 0) n0 = 0.0;
else {
t0 *= t0;
n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0);
}
float t1 = 0.6 - x1*x1 - y1*y1 - z1*z1;
if(t1 < 0) n1 = 0.0;
else {
t1 *= t1;
n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1);
}
float t2 = 0.6 - x2*x2 - y2*y2 - z2*z2;
if(t2 < 0) n2 = 0.0;
else {
t2 *= t2;
n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2);
}
float t3 = 0.6 - x3*x3 - y3*y3 - z3*z3;
if(t3 < 0) n3 = 0.0;
else {
t3 *= t3;
n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3);
}
// Add contributions from each corner to get the final noise value.
// The result is scaled to stay just inside [-1,1]
fac += (32.0*(n0 + n1 + n2 + n3))*pwr;
pwr *= pwHL;
xin *= lacunarity;
yin *= lacunarity;
zin *= lacunarity;
}
} else if ( Dim == 4 ) {
// 4D simplex noise, better simplex rank ordering method 2012-03-09
fac = 0;
float x=Pos[0]*Scale, y=Pos[1]*Scale, z=Pos[3]*Scale, w=t*Scale;
for(int p=0; p float n0, n1, n2, n3, n4; // Noise contributions from the five corners
// Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in
float s = (x + y + z + w) * F4; // Factor for 4D skewing
int i = fastfloor(x + s);
int j = fastfloor(y + s);
int k = fastfloor(z + s);
int l = fastfloor(w + s);
float t = (i + j + k + l) * G4; // Factor for 4D unskewing
float X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space
float Y0 = j - t;
float Z0 = k - t;
float W0 = l - t;
float x0 = x - X0; // The x,y,z,w distances from the cell origin
float y0 = y - Y0;
float z0 = z - Z0;
float w0 = w - W0;
// For the 4D case, the simplex is a 4D shape I won't even try to describe.
// To find out which of the 24 possible simplices we're in, we need to
// determine the magnitude ordering of x0, y0, z0 and w0.
// Six pair-wise comparisons are performed between each possible pair
// of the four coordinates, and the results are used to rank the numbers.
int rankx = 0;
int ranky = 0;
int rankz = 0;
int rankw = 0;
if(x0 > y0) rankx++; else ranky++;
if(x0 > z0) rankx++; else rankz++;
if(x0 > w0) rankx++; else rankw++;
if(y0 > z0) ranky++; else rankz++;
if(y0 > w0) ranky++; else rankw++;
if(z0 > w0) rankz++; else rankw++;
int i1, j1, k1, l1; // The integer offsets for the second simplex corner
int i2, j2, k2, l2; // The integer offsets for the third simplex corner
int i3, j3, k3, l3; // The integer offsets for the fourth simplex corner
// simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
// Many values of c will never occur, since e.g. x>y>z>w makes x < z, y < w and x < w
// impossible. Only the 24 indices which have non-zero entries make any sense.
// We use a thresholding to set the coordinates in turn from the largest magnitude.
// Rank 3 denotes the largest coordinate.
i1 = rankx >= 3 ? 1 : 0;
j1 = ranky >= 3 ? 1 : 0;
k1 = rankz >= 3 ? 1 : 0;
l1 = rankw >= 3 ? 1 : 0;
// Rank 2 denotes the second largest coordinate.
i2 = rankx >= 2 ? 1 : 0;
j2 = ranky >= 2 ? 1 : 0;
k2 = rankz >= 2 ? 1 : 0;
l2 = rankw >= 2 ? 1 : 0;
// Rank 1 denotes the second smallest coordinate.
i3 = rankx >= 1 ? 1 : 0;
j3 = ranky >= 1 ? 1 : 0;
k3 = rankz >= 1 ? 1 : 0;
l3 = rankw >= 1 ? 1 : 0;
// The fifth corner has all coordinate offsets = 1, so no need to compute that.
float x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords
float y1 = y0 - j1 + G4;
float z1 = z0 - k1 + G4;
float w1 = w0 - l1 + G4;
float x2 = x0 - i2 + 2.0*G4; // Offsets for third corner in (x,y,z,w) coords
float y2 = y0 - j2 + 2.0*G4;
float z2 = z0 - k2 + 2.0*G4;
float w2 = w0 - l2 + 2.0*G4;
float x3 = x0 - i3 + 3.0*G4; // Offsets for fourth corner in (x,y,z,w) coords
float y3 = y0 - j3 + 3.0*G4;
float z3 = z0 - k3 + 3.0*G4;
float w3 = w0 - l3 + 3.0*G4;
float x4 = x0 - 1.0 + 4.0*G4; // Offsets for last corner in (x,y,z,w) coords
float y4 = y0 - 1.0 + 4.0*G4;
float z4 = z0 - 1.0 + 4.0*G4;
float w4 = w0 - 1.0 + 4.0*G4;
// Work out the hashed gradient indices of the five simplex corners
int ii = i & 255;
int jj = j & 255;
int kk = k & 255;
int ll = l & 255;
int gi0 = perm[ii+perm[jj+perm[kk+perm[ll]]]] % 32;
int gi1 = perm[ii+i1+perm[jj+j1+perm[kk+k1+perm[ll+l1]]]] % 32;
int gi2 = perm[ii+i2+perm[jj+j2+perm[kk+k2+perm[ll+l2]]]] % 32;
int gi3 = perm[ii+i3+perm[jj+j3+perm[kk+k3+perm[ll+l3]]]] % 32;
int gi4 = perm[ii+1+perm[jj+1+perm[kk+1+perm[ll+1]]]] % 32;
// Calculate the contribution from the five corners
// MJA because I couldn't get OSL to compile an array of structs
// the 4-vectors for the gradients are split into an array with
// regular vectors (x,y,z) in grad4v and an array of floats
// (for the w component) in grad4t
float t0 = 0.6 - x0*x0 - y0*y0 - z0*z0 - w0*w0;
if(t0 < 0) n0 = 0.0;
else {
t0 *= t0;
n0 = t0 * t0 * dot(grad4v[gi0], grad4t[gi0], x0, y0, z0, w0);
}
float t1 = 0.6 - x1*x1 - y1*y1 - z1*z1 - w1*w1;
if(t1 < 0) n1 = 0.0;
else {
t1 *= t1;
n1 = t1 * t1 * dot(grad4v[gi1], grad4t[gi1], x1, y1, z1, w1);
}
float t2 = 0.6 - x2*x2 - y2*y2 - z2*z2 - w2*w2;
if(t2 < 0) n2 = 0.0;
else {
t2 *= t2;
n2 = t2 * t2 * dot(grad4v[gi2], grad4t[gi2], x2, y2, z2, w2);
}
float t3 = 0.6 - x3*x3 - y3*y3 - z3*z3 - w3*w3;
if(t3 < 0) n3 = 0.0;
else {
t3 *= t3;
n3 = t3 * t3 * dot(grad4v[gi3], grad4t[gi3], x3, y3, z3, w3);
}
float t4 = 0.6 - x4*x4 - y4*y4 - z4*z4 - w4*w4;
if(t4 < 0) n4 = 0.0;
else {
t4 *= t4;
n4 = t4 * t4 * dot(grad4v[gi4], grad4t[gi4], x4, y4, z4, w4);
}
// Sum up and scale the result to cover the range [-1,1]
fac += (27.0 * (n0 + n1 + n2 + n3 + n4))*pwr;
pwr *= pwHL;
x *= lacunarity;
y *= lacunarity;
z *= lacunarity;
w *= lacunarity;
}
}
}
Example node setup
The sample images were created with the following node setup:
