Page History

/// Oiginal notice:

/* 
 * srdnoise23, Simplex noise with rotating gradients
 * and a true analytic derivative in 2D and 3D.
 *
 * This is version 2 of srdnoise23 written in early 2008.
 * A stupid bug was corrected. Do not use earlier versions.
 *
 * Author: Stefan Gustavson, 2003-2008
 *
 * Contact: stefan.gustavson@gmail.com
 *
 * This code was GPL licensed until February 2011.
 * As the original author of this code, I hereby
 * release it into the public domain.
 * Please feel free to use it for whatever you want.
 * Credit is appreciated where appropriate, and I also
 * appreciate being told where this code finds any use,
 * but you may do as you like.
 */
 
/// Translated and modified by Ivan Mavrov, Chaos Group Ltd. 2016
/// Contact: ivan.mavrov@chaosgroup.com

/// 3D Simplex flow noise implementation.
/// Returned values are in the [0, 1] range.
float simplexFlowNoise3D(point position, float flow, output vector dNoise) {
 	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
 	};
 	
 	// Gradient component that leads to a vector of length sqrt(2).
 	// float a = sqrt(2)/sqrt(3);
 	float a = 0.81649658;
 	
 	vector gradientUBase[16] = {
 		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(-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 gradientVBase[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), 
 		vector( 0.0,  1.0, -1.0), vector( 0.0, -1.0, -1.0)
 	};
 	
 	// Helper function to compute the rotated gradient.
 	vector getGradient(int index, float sinTheta, float cosTheta) {
 		int safeIndex = index % 16;
 		vector gradientU = gradientUBase[safeIndex];
 		vector gradientV = gradientVBase[safeIndex];
 		return cosTheta * gradientU + sinTheta * gradientV;
 	}
 	
 	// Skewing factors for the 3D simplex grid.
 	// float F3 = 1.0 / 3.0;
 	// float G3 = 1.0 / 6.0;
 	float F3 = 0.333333333;
 	float G3 = 0.166666667;
 	float G3a = 2.0 * G3;
 	float G3b = 3.0 * G3 - 1.0;

 	// Skew the input space to determine the simplex cell we are in.
 	float skew = (position[0] + position[1] + position[2]) * F3;
 	point skewedPosition = position + skew;
 	point skewedCellOrigin = floor(skewedPosition);

 	// Unskew the cell origin
 	float unskew = (skewedCellOrigin[0] + skewedCellOrigin[1] + skewedCellOrigin[2]) * G3;
 	point cellOrigin = skewedCellOrigin - unskew;
 	
 	// The offset from the cell's origin a.k.a point 0
 	vector offset0 = position - cellOrigin;

 	// The second point offset from the cell origin in skewed space.
 	vector skewedOffset1;
 	// The third point offset from the cell origin in skewed space.
 	vector skewedOffset2;

 	if (offset0[0] >= offset0[1]) {
 		if (offset0[1] >= offset0[2]) {
 			// X Y Z order
 			skewedOffset1 = vector(1, 0, 0);
 			skewedOffset2 = vector(1, 1, 0);
 		} else if (offset0[0] >= offset0[2]) {
 			// X Z Y order
 			skewedOffset1 = vector(1, 0, 0);
 			skewedOffset2 = vector(1, 0, 1);
 		} else { 
 			// Z X Y order
 			skewedOffset1 = vector(0, 0, 1);
 			skewedOffset2 = vector(1, 0, 1);
 		}
 	} else {
 		if (offset0[1] < offset0[2]) {
 			// Z Y X order
 			skewedOffset1 = vector(0, 0, 1);
 			skewedOffset2 = vector(0, 1, 1);
 		} else if (offset0[0] < offset0[2]) {
 			// Y Z X order
 			skewedOffset1 = vector(0, 1, 0);
 			skewedOffset2 = vector(0, 1, 1);
 		} else {
 			// Y X Z order
 			skewedOffset1 = vector(0, 1, 0);
 			skewedOffset2 = vector(1, 1, 0);
 		}
 	}

 	// A step of (1, 0, 0) in skewed space means a step of (1 - G3, -G3, -G3) in regular space.
 	// A step of (0, 1, 0) in skewed space means a step of (-G3, 1 - G3, -G3) in regular space.
 	// A step of (0, 0, 1) in skewed space means a step of (-G3, -G3, 1 - G3) in regular space.

 	// The offset from point 1 in regular space.
 	vector offset1 = offset0 - skewedOffset1 + G3;
 	
 	// The offset from point 2 in regular space.
 	vector offset2 = offset0 - skewedOffset2 + G3a;

 	// The offset from point 3 in regular space.
 	vector offset3 = offset0 + G3b;

 	// Wrap the integer indices at 256, to avoid indexing perm[] out of bounds.
 	int i = int(skewedCellOrigin[0]) & 255;
 	int j = int(skewedCellOrigin[1]) & 255;
 	int k = int(skewedCellOrigin[2]) & 255;
 	
 	// Sine and cosine for the gradient rotation angle.
 	float sinTheta = 0.0;
 	float cosTheta = 0.0;
 	sincos(M_2PI * flow, sinTheta, cosTheta);
 	
 	// Calculate the contribution from the four points.
 	float t0 = 0.6 - dot(offset0, offset0);
 	float t02 = 0.0;
 	float t04 = 0.0;
 	vector gradient0 = vector(0.0);
 	float n0 = 0.0;
 	if (t0 < 0.0) {
 		t0 = 0.0;
 	} else {
 		t02 = t0 * t0;
 		t04 = t02 * t02;
 		gradient0 = getGradient(perm[i + perm[j + perm[k]]], sinTheta, cosTheta);
 		n0 = t04 * dot(gradient0, offset0);
 	}

 	float t1 = 0.6 - dot(offset1, offset1);
 	float t12 = 0.0;
 	float t14 = 0.0;
 	vector gradient1 = vector(0.0);
 	float n1 = 0.0;
 	if (t1 < 0.0) {
 		t1 = 0.0;
 	} else {
 		t12 = t1 * t1;
 		t14 = t12 * t12;
 		gradient1 = getGradient(perm[i + int(skewedOffset1[0]) + perm[j + int(skewedOffset1[1]) + perm[k + int(skewedOffset1[2])]]], sinTheta, cosTheta);
 		n1 = t14 * dot(gradient1, offset1);
 	}

 	float t2 = 0.6 - dot(offset2, offset2);
 	float t22 = 0.0;
 	float t24 = 0.0;
 	vector gradient2 = vector(0.0);
 	float n2 = 0.0;
 	if(t2 < 0.0) {
 		t2 = 0.0;
 	} else {
 		t22 = t2 * t2;
 		t24 = t22 * t22;
 		gradient2 = getGradient(perm[i + int(skewedOffset2[0]) + perm[j + int(skewedOffset2[1]) + perm[k + int(skewedOffset2[2])]]], sinTheta, cosTheta);
 		n2 = t24 * dot(gradient2, offset2);
 	}
 	
 	float t3 = 0.6 - dot(offset3, offset3);
 	float t32 = 0.0;
 	float t34 = 0.0;
 	vector gradient3 = vector(0.0);
 	float n3 = 0.0;
 	if(t3 < 0.0) {
 		t32 = 0.0;
 	} else {
 		t32 = t3 * t3;
 		t34 = t32 * t32;
 		gradient3 = getGradient(perm[i + 1 + perm[j + 1 + perm[k + 1]]], sinTheta, cosTheta);
 		n3 = t34 * dot(gradient3, offset3);
 	}

 	// Accumulate the contributions from each point to get the final noise value.
 	// The result is scaled to return values in the range [-1,1].
 	float noiseValue = 32.0 * (n0 + n1 + n2 + n3);

 	// Compute noise derivative.
 	//  *dnoise_dx = -8.0f * t02 * t0 * x0 * dot(gx0, gy0, gz0, x0, y0, z0) + t04 * gx0;
 	//  *dnoise_dy = -8.0f * t02 * t0 * y0 * dot(gx0, gy0, gz0, x0, y0, z0) + t04 * gy0;
 	//  *dnoise_dz = -8.0f * t02 * t0 * z0 * dot(gx0, gy0, gz0, x0, y0, z0) + t04 * gz0;
 	//  *dnoise_dx += -8.0f * t12 * t1 * x1 * dot(gx1, gy1, gz1, x1, y1, z1) + t14 * gx1;
 	//  *dnoise_dy += -8.0f * t12 * t1 * y1 * dot(gx1, gy1, gz1, x1, y1, z1) + t14 * gy1;
 	//  *dnoise_dz += -8.0f * t12 * t1 * z1 * dot(gx1, gy1, gz1, x1, y1, z1) + t14 * gz1;
 	//  *dnoise_dx += -8.0f * t22 * t2 * x2 * dot(gx2, gy2, gz2, x2, y2, z2) + t24 * gx2;
 	//  *dnoise_dy += -8.0f * t22 * t2 * y2 * dot(gx2, gy2, gz2, x2, y2, z2) + t24 * gy2;
 	//  *dnoise_dz += -8.0f * t22 * t2 * z2 * dot(gx2, gy2, gz2, x2, y2, z2) + t24 * gz2;
 	//  *dnoise_dx += -8.0f * t32 * t3 * x3 * dot(gx3, gy3, gz3, x3, y3, z3) + t34 * gx3;
 	//  *dnoise_dy += -8.0f * t32 * t3 * y3 * dot(gx3, gy3, gz3, x3, y3, z3) + t34 * gy3;
 	//  *dnoise_dz += -8.0f * t32 * t3 * z3 * dot(gx3, gy3, gz3, x3, y3, z3) + t34 * gz3;
 	
	dNoise dNoise = t02 * t0 * offset0 * dot(gradient0, offset0) + t04 * gradient0;
 	dNoise += t12 * t1 * offset1 * dot(gradient1, offset1) + t14 * gradient1;
 	dNoise += t22 * t2 * offset2 * dot(gradient2, offset2) + t24 * gradient2;
 	dNoise += t32 * t3 * offset3 * dot(gradient3, offset3) + t34 * gradient3; 	
 	dNoise *= -8.0;
 	
 	// Scale noise derivative to match the noise scaling.
 	//dNoise *= (32.0 / 2.0); 
 	
 	// Scale to the [0, 1] range.
 	return 0.5 * noiseValue + 0.5;
}

/// 3D Fractal Simplex flow noise implementation.
/// Returned values are in the [-1, 1] range.
float fractalSimplexFlowNoise3D(
 	point position,
 	float flow,
 	float lacunarity,
 	float flowRate,
 	float gain,
 	float advect,
 	int octaveCount
) {
 	float noiseValue = 0.0;
 	
 	float flowValue = flow;
 	
 	float amplitude = 1.0;
 	float advectionAmount = advect;
 	
 	for (int octave = 0; octave < octaveCount; ++octave) {
 		vector flownoise3DGradient = vector(0.0);
 		float noiseOctave = amplitude * (simplexFlowNoise3D(position, flowValue, flownoise3DGradient) - 0.5);
 		noiseValue += noiseOctave;
 		
 		position -= advectionAmount * noiseOctave * flownoise3DGradient;

 		position *= lacunarity;
 		flowValue *= flowRate;
 		amplitude *= gain;
 		advectionAmount *= advect;
 	}
 	
 	return noiseValue;
}

shader FractalSimplexFlowNoise 
 	[[ string description = "Fractal Simplex flow noise" ]] 
(
 	float start_frequency = 1.0
 		[[ string description = "Initial sampling position multiplier that affects the overall granularity." ]],
 	vector start_offset = vector(0.0)
 		[[ string description = "Offsets the initial sampling position effectively shifting the pattern in the specified direction." ]],
 		
 	float flow = 1.0
 		[[ string description = "The coordinate of a special noise dimension with a period of 1 that naturally evolves the noise to animate it instead of sliding a 3D slice throught the noise space." ]],
 	
 	float lacunarity = 2.0
 		[[ string description = "Position (frequency) multiplier per iteration." ]],
 	float flow_rate = 1.0
 		[[ string description = "Flow coordinate multiplier per iteration." ]],
 	float gain = 0.5
 		[[ string description = "Amplitude multiplier per iteration." ]],
 	float advect = 0.5
 		[[ string description = "Both initial advection amount and advection multiplier per iteration." ]],
 	int octaves = 8
 		[[ string description = "Number of fractal iterations." ]],

 	float attenuation = 1.0
 		[[ string description = "The power of the falloff applied to the final result." ]],
 	
 	output color result = 0.0
) { 	
 	point objectPosition = transform("object", P);
 	point startPosition = start_frequency * objectPosition - start_offset;

 	float noiseValue = fractalSimplexFlowNoise3D(startPosition, flow, lacunarity, flow_rate, gain, advect, octaves);
 	
 	noiseValue = 0.5 * noiseValue + 0.5;
 	noiseValue = pow(noiseValue, attenuation);
 	
 	result = color(noiseValue);
}