Version 0.0.4 of my new tree add-on is focused on defining the crown shape with groups of objects
usability enhancements
In the previous version of the add-on the markers or endpoints that represent the resources that the growing branches compete for were distributed inside an invisible ellipsoid. The resulting distribution could me made visible with the show markers option but from a usability point of view that was not ideal. likewise the su bstraction of markers from the distribution for those markers that were in the set of currently selected objects was less than intuitive, the more so since producing a tree would deselectsll objects forcing the user to select them again when trying for a different tree shape.
In this new version we allow for different way to define a crown shape using groups. If you select the use groups option you can choose three different groups:
the Crown group
any markers generated will lie inside the objects that are part of this group.
the Shadow group
any markers generated will not lie in any objects in this group. This group is optional.
the Exclusion group
growing branches will not penetrate into any objects in this group but stop. This group is also optional.
The shape of the marker cloud is determined by the combined shapes of the objects in the crown group minus those in the shadow group. If you use a shadow group it must be different from the crown group otherwise no markers would be produced at all (although there is a safeguard in the add-on to prevent this). If you use an exclusion group this should differ from the crown group as well although it might be the same as the shadow group.
workflow
The workflow now becomes pretty straight forward:
create a crowngroup, for example a bunch of different icospheres,
assign anything casting a shadow to a shadow group (here a box),
optional: assign anything impenetrable to an exclusion group (or reuse the shadow group),
position the 3d cursor at the place where the bole of the tree should start,
click update tree (here shown with the parameter settings that produce the tree in the screenshot),
Of course you can alter the other options as well. Don't forget to disable any objects for rendering that you used just to define a group otherwise you won't see much of your tree.
Code availability
This release is version 0.0.4 and avaible as a zipfile on GitHub
Blender already has a few add-ons to generate plant or tree objects, notably IvyGen and Sapling and recently I added a general L-system add-on into the mix. In this post I present a new add-on that creates trees based on the space colonization algorithm by Adam Runions.
A new kind of tree
This add-on now has a professional cousin available on Blender Market. It offers new functionality like controlling crown shapes with the grease pencil and is bundled with quite a few ready to use trees that can be used as is or function as a starting point. A YouTube playlist with examples/tutorials is available as well.
The GitHub repository now serves version 0.0.4 which makes it possible to define the shape of the crown with groups of objects. Details can be found in this follow-up article.
The GitHub repository now serves version 0.0.3 which allows more than one leaf per internode and has a control to determine how much these leaves are clumped. Note that the oft maximum of leaves is still 1.0, if you want more, click and alter the value by hand.
The GitHub repository now serves version 0.0.2 which adds the ability the enable the generation of new endpoint markers while the tree is generated. This will allow for bushier trees with extra small twigs on the main branches
[Some additional sample trees are near the end of this article] Because IvyGen is limited to climbing shrubs and Sapling requires a fair amount of user interaction to get nice results, I thought it would be interesting to explore the possibilities of generating trees with the space colonization algorithm as described by Adam Runions. This algorithm models the growth of a tree by letting branches compete for resources (think light) as they grow. The distribution and density of these resources control the final shape and the algorithm lends itself easily to adapt to environmental influences like walls blocking light etc. a feature that comes in handy when doing architectural modelling for example. In the first version of this add-on many features are already available but it is not yet a finished tool. Any remsrks and criticisms are welcome of course. This article describes how to install and use the add-on but does not give an in depth account of its design or underlying principles, that is left for a later article. I do present some example images with their parameter settings to get you started. The code is available on github. Just download the latest zip from the release directory and install the add-on with File->user preferences->addons install from file by pointing it to the zip file. You can then toggle the checkbox to use it (It is in the Add Mesh category as SCA Tree Generator). Remebr to uninstall the add-on first if you installed an earlier version by clicking the remove button.
Creating a tree
To add a tree, position your 3d cursor where you would like the tree to grow and select Add->Mesh->Add Tree to Scene. A small tree skeleton will appear and in the toolbar on the left of the 3d view will be a number of options. (Press T if the tool bar isn't visible) User interaction is copied from IvyGen so changing the values of options doesn't have an immediate effect but new values are applied when you click update tree. This is necessary because generating larger trees can take quite some time and if that was attempted on each value change the add-on would become unusably slow. An overview of all the options is available with the source, here we just highlight the more important ones.(tip: before you start experimenting it might be a good idea to check the no modifiers checkbox because the tree trunk is given substance by the skin modifier and will take many seconds on larger trees} The whole idea of the space colonization algorithm is to populate a volume (an ellipsoid in our case) with markers and let a growth point (at the the 3d cursor) start growing in the direction of the nearest markers. As a growth point (called internode or branchpoint in the add-on in a not entirely botanically correct fashion) approaches a marker close enough, this marker is removed and subsequent growth is in the direction of a slightly different set of nearest markers. A branchpoint can have any number of markers in its sphere of influence but conversely a marker has only a single closest branchpoint. Because of that branches may develop sideshoots (read more on this in the original article) This means that the number of markers and their distribution are the crucial parameters in determining the shape of the final tree,(tip if you check the show markers option in the debug panel, the markers that were generated at the start of the tree generation will be visualized as small tetrahedrons. That way you get a clear picture of their distribution. This feature is used in most images that follow.)
Some example trees and their marker distributions
Note that only the parameters that differ from their defaults are shown in the picture.
Few endpoints vs. many endpoints
Near central axis vs. near surface
More endpoints near the top
Mix and match
Even shrubs are possible as shown on the right.
adding leaves
When you click the add leaves button a new panel with options will appear with controls for the leaves. The number of leaves controls the probability that an internode (branch segment) carries a leaf. The size and size variation of these leaves can be controlled, as well as the variation in orientation of the leaves. The number of connections option controls how much 'old wood' carries leaves. If you set this to 1 only new twigs will carry leaves.
Shaping the branches
When you click on the add modifiers option, three modifiers are added to the trunk skeleton (shown in the opening image of this article):
a subsurface modifier to smooth the branches,
a skin modifier to give volume to the branches and
another subsurface modifier to smooth this skin.
Note that calculating these last two modifiers may take many seconds or even minutes if the tree is complex! The thickness of the trunk and branches and the way the thickness of forking branches is combined is controlled by the Scale and Power options respectively.
Environmental influence
A final feature of this add-on that almost comes naturally with the space colonization algorithm is the ability to factor in environmental influences. For example if one side of the growing tree receives more shadow than the other side due to a wall or building, that side will grow less. We can easily mimick this by not placing any markers in a shadowed area. To accomplish this the add-on checks whether any of the markers is generates is inside any of the selected mesh objects in the scene and discards them if they are. The workflow is then:
select one or more objects (first image below)
position the 3d cursor where you want the tree to start
tweak any options
update the tree
The result (with markers) is shown in the second image below.
Note that we did not select the actual wall, but a larger duplicate (that we will not render), because the shadow influence normally extends beyond the wall and because a new branch might still pass through this exclusion zone towards markers that are above it. By chosing a wider volume we prevent those branches to pass through the actual wall.
Limitations
In this first version there will undoubtably be any number of bugs. One of the most annoying issues is that when the kill distance is chosen too small, a developing branch may terminate in a myriad of tiny segments. This will not only take an enormous amount of time when adding the skin modifier but will also cause a strange clump of hundreds of leaves at that point. Another thing is that even when not actually generating new geometry, the user interface is sometimes slow when modifiers are selected. I have a hunch this might be related to Blender keeping track of this complicated geometry in its undo cache but I am not sure. Of course there is also plenty of room for functional improvements: it would be nice if you could select more than one starting point for example to create hedgerows etc. and also have arbitrary objects instead of square leaves, but those are nice to haves for the future.
Some more examples
The following images show the same type of tree as the opening image of this article, but at more mature stages: The masked leaf texture as well as the concrete and brick textures are from cgtextures.com, the latter two converted to various types of specular and bump maps with the free shadermap commandline tool. The hdri-background is Alex Hart's 'Arboretrum in Bloom'. The most mature tree was generated with the settings shown below:
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 attempt to create an X-ray shader that highlights object edges we profit form the built-in facilities of OSL to compute derivatives.
The X-ray effect in the picture is mostly achieved by letting the meshes emit some light but more from their contours as seen from the camera.
(note that with all that transparency we need an awfull lot of samples to get a noise free result: even with 500 samples there is still noise visible in the image)
In order to determine what the contours are, we make use of OSLs built-in functions Dx() and Dy() that compute the derivatives of a function. The idea is that the derivatives for the current shading possition as seen from the camera change fastest at the contours of objects. This of course will only work reasonably well for curved objects. In the picture below we have a simple diffuse that is red when the sum of the derivatives is large. The small cube on the left has sharp edges that case an abrupt change that we cannot capture this way. The cube on the right has a bevel and a subsurface modifier added and does show edges.
shader der( point Pos = P, output vector dx =0, output vector dy =0, output float Lx = 0, output float Ly = 0, output float R = 0 ){ dx = Dx(Pos); dy = Dy(Pos); Lx= length(dx); Ly= length(dy); R = Lx + Ly; }
Sample node setup
The x-ray image at the beginning of this post was created with the following node setup:
Beside a mix of shaders the most important part is the texture coordinate that we take as input: Camera space is selected here. The add and multiplication nodes are there just to give us some control over the noise we want to mix in.
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:
Simplex noise is a simple and fast alternative to Perlin noise that scales easily to arbitray dimensions. In this post I present a an implementation for 2, 3 and 4 dimensions
Simplex Noise
as of 22/4/2013 osl supports simplex noise natively, albeit for dimensions up to 4 only and without advection. This means that the code presented here and other articles on this blog stays relevant even when Blender incorporates this new OSL version but by then it probably would make sense to use the built-in simplex noise for 3 and 4 dimensions.
Simplex Noise is a simpler and faster alternative the the ubiquitous Perlin noise (Blender calls it simply 'noise'in OSL). It was designed by Ken Perlin a decade ago and a nice article by Stefan Gustavson on the details and possible implementation is a good starter if you want to know more. In this post I present a port to OSL that is a more or less straight reimplementation of his java code.
The sampler above shows (from left to right) 2, 3 and 4D noise on a flat square and a sphere. The animated gif shows a short sequence of 4D noise (click to play from my site, Blogger doesn't show animated gifs).
The code is quite long, but well commented. I retained most of the comments from Stefans Java code as they explain rather well what is going on. To port the code pretty much all I needed to do was replace doubles by floats and to write a shader with some functions (OSL isn't object oriented so I had to seperate those out). The only difficulty was that the OSL compiler refuses to compile arrays of structs (I will investigate that further since it should be possible). Therefor, instead of defining an array of structs to represent a four dimensional gradients I split those arrays into an array of vectors and an array of floats. See the comments in the code where that is relevant. I realize that this code is quite long and especially the first part is hard to read due to the wide permutation tables so I advise to download it (click view plain) and read it in your favorite editor. (If you use Notepad++ or Ultraedit, setting syntax highlighting to javascript or C works quite well, Blenders internal editor has native highlighting support for OSL if you give the text buffer an .osl extension but I doesn't do a nice job on this example, at least not in my revision (r53600)).
/*
* 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.
*
*/
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 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;
if(Dim == 2){
// 2D simplex noise
float xin=Pos[0]*Scale, yin=Pos[1]*Scale;
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);
} else if(Dim == 3){
// 3D simplex noise
float xin=Pos[0]*Scale, yin=Pos[1]*Scale, zin=Pos[2]*Scale;
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);
} else if ( Dim == 4 ) {
// 4D simplex noise, better simplex rank ordering method 2012-03-09
float x=Pos[0]*Scale, y=Pos[1]*Scale, z=Pos[3]*Scale, w=t*Scale;
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);
}
}
In this post we tabulate some of the parameters of Gabor noise as provided in Cycles'Open Shading Language.
In a previous post I presented a small OSL node that gives access to the fairly new Gabor noise that is available in OSL. I already showed some examples of how to use it and made an attempt to implement a fabric shader based on this noise. To help visualize the parameters of this noise I made the overview shown below (click to enlarge):
Gabor noise has tree modes:
Isotropic
Anisotropic
Mixed
These modes are triggered by the values 0, 1 and 2 of the anisotropic parameter. In isotropic mode the direction parameter is ignored but in the other modes the length and direction of this vector determine the frequency and the direction of the anistropy.
The impulses parameter determines the number of features (bumps) present per unit area, whereas the bandwidth parameter determines how much detail of these bumps is visible.
In the image above the colored squares show a systematic variation of these parameters. Gabor noise has a value between [-1, 1] centered around zero (the green in the colored examples). The black and white examples show some specific patterning that is possible if you choose suitable cutoffs for the noise values. They we created using the following node setup:
Note that the add/divide nodes on the left to nothing more than map the [-1, 1] range to [0, 1]. The less than/ more than nodes in combination with the multiply node pass on a value if its between the two limits or zero otherwise.
A combination of two perpendiclar anisotropic Gabor textures is an easy way to create textures that mimic fabrics. In this post we use this to create a simple denim fabric.
In a previous post I presented a simple node to provide Gabor noise. Anisotropic Gabor noise with a large direction vector has slightly undulating parallel lines. Two of these perpendicular to each other are an easy way to create a fabric like structure.
(Click to enlarge if you cannot see the individual threads of the fabric)
The material that is presented here is probably not as versatile or fast as the fabric node I presented earlier, but it does give a somewhat more natural result (at least in my opinion) and is foremost an example of how versatile Gabor noise is. Also, because a Gabor noise node is farly high on the todo list for Cycles this means in the future we might not be dependent on custom nodes as presented here.
Example node setup
The node setup here isn't exactly pretty (click to enlarge) but the main idea is the combination of two Gabor noise nodes (show in a previous post):
As indicated in the image both nodes have a fairly large direction vector which gives small lines. Both vectors are prependicular. The resulting noise is a value that may be less than zero and adding a value to the result as we do here allows for control over the apparent gap between the threads of the fabric. A hight value gives less visible space between the threads. Note that the scale value on the left controls the scaling of the UV-space. It has an effect on the variation of the noise but not on the number of threads per centimeter, that is completely controlled byt the direction vectors. The bandwidth value controls the quality: changing it gives less or more of a 'fluffy' appearence to the threads.
